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ABSTRACT 


Equations  are  derived  for  the  species  concentrations  in  coupled 
chemical  equilibria  in  the  gas  phase.  This  is  done  both  for  the  method 
of  the  direct  minimization  of  the  free  energy  as  well  as  for  the  equilib¬ 
rium  constant  method.  The  relations  developed  are  such  as  to  allow  for 
the  inclusion  of  detailed  real  gas  effects.  The  emphasis  is  on  practical 
problems  encountered  in  the  actual  calculation  of  the  species  and  the 
thermodynamic  properties.  Expressions  are  derived  for  the  direct  cal¬ 
culation  of  the  concentration  derivatives  with  respect  to  temperature 
and  density  required  for  the  calculation  of  the  derived  properties  (i.  e. 
specific  heats,  sound  velocities,  etc.  ).  Considerable  space  is  given  to 
discussion  of  the  non-linear  numerical  methods  available  for  the  solu¬ 
tion  of  the  non-linear  equations  for  the  species  concentrations. 
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1.  INTRODUCTION 


In  this  report  we  shall  discuss  the  determination  of  the  properties  of 
a  gaseous  mixture  capable  of  changing  composition  through  chemical  reaction. 

Our  approach  will  differ  somewhat  from  the  usual  textbook  approach  in  the 
sense  that  we  shall  go  beyond  the  development  of  the  fundamental  thermodynamic 
relations  to  discuss  some  of  the  practical  problems  which  arise  in  the  actual  ^ 
calculation  of  the  thermodynamic  properties  and  compositions  for  such  mixtures.  “ 
We  shall  restrict  our  discussion  to  homogeneous  gas  phase  reactions,  the 
interested  reader  being  referred  to  the  appropriate  literature  for  extensions 
which  include  condensed  phases.® 

In  the  derivations  that  follow,  the  usual  assumption  is  made  that 
chemical  reactions  can  be  frozen  at  any  point  in  their  approach  to  equilibrium. 

In  other  words,  our  system  can  be  in  thermal  and  mechanical  but  not  chemical 
equilibrium.  The  reasonableness  of  treating  such  systems  within  the  framework 
of  thermodynamics  is  based  essentially  on  the  fact  that  the  functions  which 
describe  the  properties  of  the  system  when  it  is  not  in  chemical  equilibrium 
reduce  to  those  for  complete  thermodynamic  equilibrium  on  substitution  of  the 
equilibrium  compositions,  no  other  changes  being  required.® 

i 

The  analysis  of  the  chemical  equilibrium  problem  ultimately  results  in  a 
set  of  non-linear  equations  which  have  to  be  solved.  Since  these  equations  can¬ 
not  be  solved  in  closed  form,  numerical  procedures  for  their  solution  must  be 
developed.  These  numerical  procedures  then  form  an  integral  part  of  the  overall 
solution  of  the  problem.  For  this  reason,  a  portion  of  this  chapter  has  been 
devoted  to  the  discussion  of  appropriate  numerical  methods. 

/ 

Of  primary  importance  in  chemical  equilibrium  calculations  is  the  selec¬ 
tion  of  a  subset  of  species  concentrations  to  be  a  basis  set  (in  the  sense  of 
vector  analysis)^  for  the  mathematical  description  of  the  problem.  'There  will 
be  a  minimum  number  of  concentrations  which  span  the  space  of  all  species 
concentrations.  The  choice  of  particular  species  concentrations  for  this  basis 
set  is  somewhat  arbitrary.  The  species  whose  concentrations  have  been  selected 
to  be  that  set  will  be  referred  to  ag  reference  species,  all  other  species 
being  designated  derived  species ■  Since  the  number  of  reference  species  is 
Intimately  connected  to  the  number  of  distinct  atomic  species,  a  certain  sim¬ 
plicity  is  obtained  if  the  atoms  (and  free  electrons)  are  chosen  as  reference 
species.  This  particular  choice  of  reference  species  is  not  always  convenient 
from  a  computational  point  of  view  however.  That  problem  is  discussed  below  and 
simple  matrix  methods  are  presented  for  the  transformation  of  species  from 
derived  to  reference  species  and  vice  versa.  In  most  of  our  discussion,  we 
shall  consider  the  gaseous  mixture  to  consist  of  l  chemical  species  made  up  of 
free  electrons,  c-1  atomic  species,  and  m(=&-c)  other  species,  the  atomic  species 
together  with  the  electron  being  designated  as  reference  species  and  all  other 
constituents  being  designated  derived  species. 
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The  chemical  formula  of  any  given  species  can  be  written,  symbolically, 
as 


n  (PJ 

j-1  j  ViJ 


i  **  1 , . . « ,  1 


(1) 


where  the  P  refer  to  the  formulae  of  the  reference  species  only  (since  we  take 
these  to  be^the  atoms  and  free  electrons).  It  should  be  noted  that  the 
reference  species  have  also  been  included  among  the  .  For  these  species, 


however,  (1)  becomes  an  identity,  since  for  them  v, ,  “  5. . ,  the  Kronecker 
delta.  The  subscripts  v  represent  the  amount  of  Reference  species  j  con¬ 
tained  in  one  molecule  or J species  i. 


A  chemical  reaction  connecting  the  four  species  S^,  S^,  and  can 
be  written  J  ^ 


Vl  +  °2S2  "  a3S3  +  %S4 


(2) 


According  to  convention,  the  left  hand  members,  S..  and  S_  are  called  reactants 
while  the  right  hand  members  ,  and  S^,  are  called  products.  In  what  follows, 
it  will  be  convenient  to  replace  the  stoichiometric  coefficients  by 
coefficients  t.  whose  magnitudes  are  equal  to  those  of  the  corresponding  a., 
but  which  are  negative  for  reactants  and  positive  for  products.  Equation  *2) 
becomes 


l  t i  Si  =  0  C3) 

i«=l 


The  coefficients  and  hence  the  are  chosen  so  that  (3)  is  balanced  for 
every  reference  species.  IhiB  balancing  can  be  simply  expressed 


i“l 


V 


ij 


(4) 
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where  j  runs  over  all  the  reference  species. 

Equation  (3)  has  been  written  for  the  single  reaction  (2)  involving  the 
four  species  S^.  Its  generalization  to  simultaneous  chemical  reactions  for 
arbitrary  numbers  of  species  is  simply 


K*3!  -  0 


C5> 


where  k  runs  over  all  the  reactions  and  i  over  all  species.  Equation  (4)  is 
essentially  unchanged  by  the  addition  of  the  second  subscript  to  form  t  ,  .  It 
will  be  noted  that  t^  is  zero  when  the  ic^  reaction  does  not  contain  trie 

species.  In  other  words,  the  matrix  T  (whose  elements  are  the  t^)  i9  defined 
such  that  each  column  19  associated  with  one  species  and  each  row  with  one 
reaction  ordered  in  some  arbitrary  fashion.  Each  row  column  position  of  such 
a  matrix  is  then  associated  with  a  particular  chemical  reaction,  and  a  particu¬ 
lar  chemical  species  which  takes  part  in  that  reaction. 

I 

The  specification  of  the  thermodynamic  state  of  a  heterogeneous  mixture 
requires  the  specification  of  two  thermodynamic  variables,  e.g„,  temperature 
and  pressure,  along  with  the  specification  of  the  composition  of  each  phase. 

Since  we  shall  restrict  ourselves  to  the  gas  phase,  it  is  only  necessary  to 
specify  the  pair  of  thermodynamic  variables  along  with  the  composition  for  the 
single  phase.  The  specification  of  the  path  by  which  this  thermodynamic  state 
was  obtained  is  totally  irrelevant  to  the  specification  of  the  state.  For 
systems  capable  of  undergoing  composition  change  through  chemical  reactions, 
this  Is  equivalent  to  stating  that  the  particular  chemical  reactions  by  means  of 
which  the  equilibrium  composition  was  attained  need  not  be  stated.  In  fact, 
in  describing  any  particular  such  state,  one  is  free  to  choose  any  convenient 
set  of  possible  reactions  which  includes  all  species  of  interest.  The  ifesults 
of  all  discussions  are  then  independent  of  such  a  choice.  In  the  initial  parts 
of  what  follows  we  shall  consider  all  molecules  and  ions,  which  we  take  as 
the  derived  species,  to  be  built  up  from  their  constituent  atoms  and  free 
electrons,  which  we  take  as  the  reference  species.  This  has  computational 
advantages  particularly  for  systems  consisting  of  molecules  containing  small 
numbers  of  atoms.  Practical  considerations  related  to  the  numerical  finite— 
ness  of  computers  will  require  the  modification  of  this  approach  when  the 
concentrations  of  the  atomic  species  and  electrons  become  extremely  small.  In 
any  event,  we  shall  always  assume  that  each  chemical  reaction  equation  contains 
only  one  derived  species,  all  other  species  in  the  equation  being  reference 
species  which  in  our  discussion  we  shall  consider  to  be  atoms  and  electrons. 

The  formation  of  all  derived  species  from  the  reference  species  reyuiiB 
in  a  simplification  in  the  T  matrix.  In  particular,  a  column  which  refers  to  a 
derived  species  will  contain  zeros  except  for  that  one  row  position  corresponding 
to  the  chemical  reaction  for  the  formation  of  that  particular  derived  species, 
since  that  derived  species  cannot  appear  in  any  other  reaction  equation.  Further¬ 
more,  this  single  non-zero  matrix  position  will  contain  the  number  plus  one.  In 
what  follows  we  shall,  therefore,  redefine  the  T  matrix  to  omit  such  columns,  that 
is,  such  that  only  columns  associated  with  the  reference  species  are  included,  it 
being  understood  that  t  =  +1  is  associated  with  the  derived  species  itself.  There 
will  be  no  confusion  on  this  latter  point  since  we  shall  shortly  explicity  insert 
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t  =  +1  for  the  derived  species  in  all  the  working  relationships.  An  example 
orJsuch  a  T  matrix,  as  redefined,  appears  in  Appendix  A. 

Any  extensive  property  of  a  mixture  can  be  written  as  a  weighted  sum 
over  the  corresponding  partial  molal  quantities  for  each  species,  each  such 
quantity  being  weighted  by  the  number  of  moles  of  the  species.  The  Gibbs  free 
energy  is  such  an  extensive  property  and  since  the  partial  molal  Gibbs  free 
energy  is  just  the  chemical  potential,  the  appropriate  weighted  sum  takes  the 
form 


G 


I 


i-1 


Vi 


(6) 


y  being  the  chemical  potential  and  n^  the  number  of  moles  of  the  ifc  species. 
At  equilibrium,  the  Gibbs  free  energy  must  be  a  minimum  with  respect  to  all 
virtual  variations  consistent  with  any  constraints  on  the  system.  In  the  ab¬ 
sence  of  nuclear  transformations,  it  is  clear  that  any  variation  carried  out 
must  be  such  as  to  preserve  the  total  number  of  reference  species,  whether 

bound  or  free.  ThiB  is  a  constraint  and  leads  to  conservation  equations 
which  can  be  written  in  the  form 

l  w  -  0  J  - 1 . c  ■ f7) 

i=l 


reference  species  other  than  electrons, x.  is  the  concentration  of  the 
j  reference  species,  whether  bound  or  free,  in  units  consistent  with  the 
n^.  For  the  electrons,  the  analagous  conservation  equation  is  most  conveniently 
expressed  in  terms  of  net  charge  conservation,  i.e.x,,  in  this  case,  represents 
the  net  overall  charge  of  the  gas.  In  particular,  for  a  neutral  gas.Xj^O. 


The  problem  with  which  we  shall  be  occupied  can  thus  be  stated  as  being  the 
determination  of  the  composition  variables  n. ,  such  that  the  free  energy  given 
by  (6)  is  a  minimum  under  all  virtual  variations,  subject  to  the  constraints  (7). 
There  have  developed  two  numerical  approaches  to  the  solution  of  this  problem. 

In  one  of  these,  the  problem  is  numerically  attacked  directly  as  stated  and,  in 
fact,  composition  variables  are  sought  which  result  in  a  stationary  value  of  G, 
subject  to  the  stated  constraints.  This  method  is  referred  to  as  the  direct 
minimization  of  the  free  energy.  In  the  second  approach,  the  formalism  is 
allowed  to  proceed  further  before  the  numerical  attack  is  mounted.  Thus,  formal 
variation  of  (6)  is  carried  out,  subject  to  the  constraints  (7).  There  results 
a  series  of  equations  each  of  which  connects  composition  variables  with  the 
composition  independent  parts  of  the  chemical  potentials  for  those  species  appearing 
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in  the  equation.  Since  these  composition  independent  terms  can  be  combined 
to  define  an  equilibrium  constant,  this  second  method  is  referred  to  as  the 
equilibrium  constant  method.  We  shall  now  proceed  to  derive  the  working 
equations  for  the  two  methods.  In  what  follows,  we  shall  formally  carry  along 
in  the  derivation  terms  containing  various  departures  from  the  ideal  gas. 

Having  derived  the  required  relations,  we  shall  then  specialize  them  to  the 
ideal  gas  in  discussing  methods  for  solution.  The  detailed  way  in  which  the  non¬ 
idealities  are  calculated  is  contained  in  Appendix  B.  Modifications  in  the 
methods  of  solution  as  required  by  these  non- idealities  will  be  indicated  but 
not  discussed  in  detail. 

2.  DIRECT  MINIMIZATION  OF  THE  FREE  ENERGY3’5’6,11 

In  principle,  the  equations  associated  with  this  method  have  been  derived, 
namely,  (6)  and  (7).  For  an  ideal  gas,  these  can  be  written  more  explicitly 
in  terms  of  the  mole  fractions  of  the  species.  Thus,  the  chemical  potential 
for  a  constituent  in  a  mixture  can  be  written 


VJi(T,P) 

RT 


U*(T,P) 

RT 


lnx^  +  lti^i 


(8) 


where  x  is  the  mole  fraction  of  the  species  in  the  mixture,  (T,  P)  the 
chemical  potential  of  the  pure  species  i  at  the  Bame  temperature  and  pressure, 
and  y  the  activity  coefficient  of  species  i  in  the  solution,  y^  contains 
all  departures  from  ideality.  For  a  real  gas,  y  depends  on  T,P,  and  the 
concentrations  of  the  various  constituents.  For  an  ideal  gas  y^=  1,  in  which 
case  the  concentration  dependence  reduces  simply  to  the  natural  logarithm  of  the 
mole  fraction,  can,  of  course,  be  written 


U*(T,P) 

RT 


o  ,  „ 
UL  (T) 

RT 


InP 


(9) 


in  which  the  standard  state  is  taken  as  1  atm.  (101325. ON/m  )  and  where  p.  (T) 
depends  only  on  the  temperature.  In  what  follows,  it  will  be  convenient  to  write 
the  free  energy,  G,  in  dimensionless  form  by  dividing  it  by  RT.  Substitution 
of  (8)  into  the  dimensionless  form  of  (6)  then  yields 


G_ 

RT 


n. 


111  < 


+  InY , 


(10) 
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The  chemical  equilibrium  problem  is  then  solved,  by  the  method  of  the  direct 
minimization  of  the  free  energy,  when  the  set  of  composition  variable  n^  is 
found  which  minimizes  (10)  subject  to  the  c  equations  (7)  being  satisfied.  It 
should  be  noted  that  when  real  gas  effects  are  included,  the  y.  can  depend  on 
the  concentrations  of  all  species.  This  will  complicate  certain  of  the 
numeric^  methods  used  in  the  solution  for  the  direct  minimization  of  the  free 
energy . 


3.  THE  EQUILIBRIUM  CONSTANT  METHOD 


As  we  have  already  stated,  this  method  starts  with  the  formal  variation  of 
(6)  subject  to  (7).  The  fact  that  the  n.  appear  in  both  (6)  and  (7)  indicates 
that  the  formal  variation  of  the  n^  in  (6)  cannot  be  carried  out  independently. 
These  variations  can  be  made  independent  by  the  elimination  of  the  proper  number 
of  variables.  This  can  be  accomplished  by  the  method  of  Lagrange  multipliers^ 
in  which  a  sum  over  the  c  equations  (7) ,  each  multiplied  by  an  unknown  multiplier 


Xj ,  is  added  to  (6).  When  this  is  done  and  the  variation  carried  out,  there 
results 

6G  = 

i-1  ‘  j  =1 


£  c 

I  (tX,  +  l  ^.v.  ■>  Sn 


+ 


l 


j=l 


£ 

(I 

i*=l 


VUai 


(ID 


where  the  unknown  Lagrange  multiplies  X.  are  defined  so  as  to  eliminate  the 
proper  number  of  coefficients  of  the  6n^.  All  remaining  variations  then  become 
independent.  The  two  sets  of  variations  in  (11)  can  then  be  carried  out 
independently.  This  leads  to  the  set  of  equations 


■i  +  *  Yu  " 0 

3=1 


and 


£ 

Y 

i«l 


ij°j  “  Xj  " 


i  =  !,...,£ 


j  -  1. 


,c 


The  first  of  these  can  be  used  to  evaluate  the  Lagrangian  multiplies  X 


j 


and 
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to  determine  their  meaning. 


Since 


"EVij* 

J 


it  follows  from  (6)  that 


l  Vi 


-  -  l 

i,  j 


A  JVi  D1 


.n , 


But,  from  the  second  set  of  relations  above 


l 

i 


Vi 


X  . 
1 


so  that 


G  = 


-  I  A-X. 

•  J  J 
J 


and  -A  can  be  considered  to  be  the  contribution  of  reference  species  j 
(whether  bound  or  free)  to  the  free  energy.  Since 


v  ax  / 


it  follows  that  -A.  is,  in  effect,  the  chemical  potential  of  the  reference 
species  whether  botlnd  or  free. 

Equation  (11)  contains  the  Z  variations  Sn..  The  fact  that  these  are 
subject  to  the  c  constraints  (7),  means  that  there  must  exist  Jt-c  variables 
which  can  be  varied  independently.  This  is  equivalent  to  the  statement  that 
there  exists  at  least  one  get  of  2-c  chemical  reactions  for  the  attainment  of 
the  thermodynamic  state.  One  such  set  of  reactions  can  be  obtained  by  considering 
each  derived  species  to  be  a  product  and  requiring  it  to  be  built  up  from  its 
constituent  reference  species,  the  latter  considered  to  be  reactants.  One  ob¬ 
tains,  thereby,  Z-c  chemical  reactions  for  the  derived  species  (along  with  c 
superflous  identity  chemical  reactions  for  the  reference  species).  Clearly 
any  other  set  of  reactions  containing  more  than  Z-c  reactions  must  contain 
redundancies . 

Having  reduced  the  number  of  reactions  to  Z-c,  one  can  obtain  an  equal 


7 


AEOC-TR-71  -52 


number  of  variables  by  defining  a  variable  to  go  with  each  reaction.  This  is 
the  sense  of  de  Donder’s1^  introduction  of  the  degree  of  reaction  variable. 
This  variable  is  defined  in  terms  of  the  change  in  the  number  of  moles  of 
a  species  produced  by  a  particular  chemical  reaction.  Thus,  the  change  in  the 
number  of  moles  of  species  j  as  a  result  of  the  ifc  reaction  is  given  by 

i  •  i  i> 


where  £  is  the  degree  of  reaction  of  chemical  reacting  i.  By  definition,  the 
same  value  of  must  apply  to  each  species  in  the  i  reaction,  the  scale 
for  each  such  |gecies  being  given  by  t,^.  The  total  change  In  the  number  of 
moles  of  the  j  species  for  the  mixture  can  be  calculated  by  summing  these 
changes  over  all  reactions.  This  yieldB 


V5i 


(12) 


which  can  be  used  in  (11)  to  demonstrate  the  Independence  of  the  variations 
Thus,  substitution  of  (12)  in  (11)  yields 


6G 


£ 

l 


i"l 


k 


+ 


c 

l 

J-l 


£ 

(I 


l  vij  ni 

i-1  J 


-  V 


6A, 


0 


I  Uifcik6ek  +  l  A  V  t  fie 
i.k  i,j ,k  2  -k  k 


c  £ 

+  I  Cl  ui1n1  “  x.i) 

j=l  i«=l  22  2  2 


But,  according  to  the  generalized  form  of  (4), 
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I 


,0 


so  that 


SG 


1Jk'itik*tk  +  U  -  V 


6A 


0 


The  fact  that  the  A  have  disappeared  from  the  first  term  on  the  right  means 
that  the  variations^have  been  uncoupled.  Independent  variation  of  the  terms 
in  this  last  equation  then  leads  to  the  relation 


j  ~  !»*■•»  i~c 


(13) 


(the  so-called  equations  of  reaction  equilibrium)  as  well  as  the  c  equations 
(7).  The  Z- c  equations  (13)  along  with  the  c  relations  (7)  also  completely 
specify  the  system.  Equations  (13)  will  form  the  basis  for  the  derivation 
of  the  working  equations  of  the  equilibrium  constant  method. 

Let  us  again  consider  the  derived  species  tD  be  formed  only  from 
reference  species.  As  previously  mentioned,  for  such  chemical  reactions  the 
coefficient  of  the  derived  species  must  be  plus  unity  so  that  (13)  can  be 
written 


t: 

+  Mil 
1=1  j  J 


i  =  1, . . . ,i-c 


(14) 


where  i  now  refers  only  to  derived  species  and  j  referee  only  to  reference 
species  and  where  there  no  longer  remains  the  possibility  for  confusion  with 
regard  to  the  removal  from  the  T  matrix  of  the  columns  corresponding  to  the 
derived  species. 

Equations  (14)  are  valid  for  all  values  of  temperature  and  pressure 
for  which  the  system  is  in  thermodynamic  equilibrium.  The  explicit  dependence 
of  these  equations  on  composition  for  the  ideal  gas  becomes  clear  on  the 
substitution  of  (9)  into  (8),  setting  y  .  =  1,  and  substituting  the  result 
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into  (14) .  We  choose,  however,  to  carry  along  the  activity-  coefficient 
for  the  present,  (14)  then  becomes 


o 


RT 


+ 


+  InP  +  It. .  InP  +  lnx  +  It. , lnx  +  lny  +  f 

*  i  j  iJ  j  i  i  y  J 


0  (15) 


The  first  two  terms  depend  only  on  the  temperature.  It  i|^customary  to  com¬ 
bine  such  terms  into  311  equilibrium  constant  for  the  I  reaction  (and 
hence  for  the  iC  derived  species)  by  the  definition 


-  InK. 
i 


Lx  +  It 

RT  j  ij  RT 


(16) 


Substitution  Into  (15)  leads  to 


-  InK, 


lnP(Ztij4+ixi+  l  lnx  1J-flnr1+  £  In 


which  can  be  written 


Ki  = 

1 

lj.  a 

- i 
u. 

ft 

t_u 

_ 1 

•  r  -*• 

1 - 

•n 

u 

X 

t=;  •<-* 

i _ 

L  J 

L  J 

(17) 


where  a>.  =  [It.  ,+l]  is  the  Increase  In  the  number  of  particles  In  the 
reaction  in  going  from  derived  to  reference  species.  In  what  follows  it 
will  be  convenient  to  introduce  an  activity  coefficient  for  the  reaction 
by  defining 
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(17)  can  then  be  written 


Vi 


p  i  nx. 

b  j 


(17  a) 


We  shall  new  restrict  ourselves  to  the  ideal  gas,  taking  y.  ■  1  as  required. 
Real  gas  effects  can  be  included  through  the  introduction  of  explicit 
expressions  for  y^  as  described  in  the  appendix. 

For  the  ideal  gas,  (17)  can  be  rewritten  so  that  the  mole  fraction  of 
the  derived  species  appears  on  the  left  and  only  reference  species  appear  on 
the  right.  Thus,  by  transposition  (17)  becomes,  for  y^  *»  1 


x  ^K-P^inxj  1J  (17b) 

•i  i  j  '• 


Since,  for  reactants  t. .  ^v..,  (17b)  can  now  be  simply  written  in  terms  of 
the  stoichiometric  coerf iciente ,  i.e. 


"i  Vij 
x,  -  k  p  nXl  J 

11  j  3 


while  the  activity  coefficient  for  the  reaction  becomes 


(18) 


h 


i 


Equation  (18)  will  form  the  basis  for  computations  within  the  framework 
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of  the  equilibrium  constant  method.  Because  we  have  chosen  to  write  the 
chemical  reactions  entirely  in  terms  of  atoms  and  the  electron  as  reference 
species ,  (18)  contains  v.  in  an  unambiguous  manner  with  regard  to  sign.  In 
what  follows,  we  shall  also  continue  to  use  t  on  occasion  in  order  to 
maintain  complete  generality.  1 


Equation  (18)  is  appropriate  for  calculations  in  which  pressure  is  a 
thermodynamic  state  variable.  If,  for  some  reason,  partial  pressures  are 
preferred  over  mole  fractions,  use  can  be  made  of  the  definition  P  =x.P  where 
P  is  the  partial  pressure  and  (18)  is  then  written  3 


v 


Ki  n.  pj 

3  3 


ij 


Quite  often  density  is  a  more  convenient  state  variable  than  pressure. 
Conversion  of  (18)  to  the  corresponding  density  form  requires  the  explicit 
use  of  the  equation  of  state  which  connects  pressure  and  density.  Thus 
since  P  is  expressed  in  bars  (or  atmospheres) ,  this  can  be  written 

X  -  k.(p/po)  1  lix  .  *3 

j  J 


where  P  is  understood  to  refer  to  the  pressure  at  specified  reference 

conditions  (usually  T*  273.15  K  and  P«  1  bar  (or  1  atm.)).  The  equation  of 

state  can  be  written  P  =  pZRT  where  Z  iB  the  compressibility  factor.  For 

the  one  component  ideal  gas,  Z  is  unity  with  departures  from  unity  being 

due  entirely  to  non- ideality.  We  shall  use  the  same  form  for  the  equation 

of  state  for  the  reacting  multicomponent  gas.  In  that  case,  it  Is 

convenient  to  let  1/p  be  the  volume  per  mole  of  the  reaction  mixture  so 

that,  for  the  ideal  gas,  Z  becomes  the  number  of  moles,  a  quantity  which 

depends  on  the  thermodynamic  state  parameters.  Naturally,  for  the  multi- 

component  real  gas,  Z  also  contains  the  effect  of  non-ideality.  New,  the 

equations  of  state  are  P  *  pZRT  for  the  conditions  of  interest  and  F  "p  Z  RT 
i  *  o  o  o  o 

at  the  reference  state  so  that 


P 

P 


T 

T 


(18)  can  then  be  written 
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where  C  =  x.  Z/Z  is  the  concentration  of  species  j  in  moles  per  mole  of 
equilibrated  gas°at  the  reference  conditions.  It  should  be  noted  that  the 
reference  state  is  being  used  for  purposes  of  scaling  all  results  and  not 
as  a  standard  state  in  the  usual  chemical  sense.  The  use  of  this  particular 
reference  state  may  be  inconvenient  in  some  systems  in  which  case,  if 
solutions  are  still  desired  at  specified  densities,  either  other  definitions 
of  the  variables  subscripted  zero  should  be  used  or  calculations  at  the  desired 
densities  carried  out  in  a  pressure  formulation  for  which  no  such  reference 
state  is  needed  with  an  iteration  inserted  so  that  the  particular  pressure 
which  corresponds  to  the  desired  density  is  determined. 


Since,  in  what  follows  a  reference  state  will  always  be  used  for 
scaling,  the  dropping  of  the  Z  should  result  in  no  confusion.  With  this  in 
mind,  Z/Z  can  be  replaced  by  2  and  the  concentrations  defined  C^=x^Z,  it 
being  understood  that  these  are  with  respect  to  one  mole  in  the  reference 
state. 


It  is  convenient  to  combine  the  factors  which  depend  only  on  temp¬ 
erature  by  defining  an  equilibrium  constant 


K. 

x 


K. 

i 


( 


T_  u> 

T  } 

o 


i 


One  obtains,  thereby. 


Ci  “  Ki(p/po)  1 

X  ^  A  J 


(19) 


For  consistency  (7)  should  be  expressed  in  the  same  units  as  (19).  Thus 
(7)  becomes 


I  vij  Cl  -  Xj  -  0  3-1 . .  (20) 

i*l 
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where  Xj  is  now  moles  per  mole  of  gas  at  the  reference  state.  Since 
for  the-*  electrons ,  C20)  represents  the  net  charge  conservation  equation 
(with  positive  and  negative  charges  balancing  each  other),  for  them  v 
(20)  is  replaced  by  t^.. 


I 


in 


Equations  (19)  and  (20)  are  the  working  equations  for  the  equilibrium 
constant  method.  These  equations  consititute  a  set  of  £  equations  for  the 
concentrations  of  all  of  the  species,  (i.e.  both  the  reference  and  derived 
species).  As  mentioned  earlier,  these  £  concentrations  are  completely 
determined  once  the  concentrations  of  the  c  reference  species  have  been 
specified.  The  problem  can  be  viewed,  therefore,  as  requiring  the  solution 
of  the  c  simultaneous  equations  (20)  subject  to  the  £-c  conditions  (19).  The 
meaning  of  this  view  of  the  problem  will  become  somewhat  more  transparent 
when  expressions  are  derived  below  for  the  calculation  of  the  derivatives  of 
the  concentrations.  The  validity  of  this  view  can  be  demonstrated  on 
substituting  (19)  into  (20).  This  yields  the  set  of  equations 


J  "  1, • • . , c 


(21) 


which  is  a  set  of  c  equations  in  c  unknowns,  i.e.  the  c  reference  species. 

Equations  (21)  are  a  set  of  highly  nonlinear  equations.  The  nonlinearity 
arises  both  from  the  exponents  v. .  on  the  unknown  variables  as  well  as  from 
the  appearance  of  products  of  the'1  unknowns.  The  question  of  the  existence  of 
multiple  solutions  naturally  arises  in  such  cases.  Such  problems  are  beyond 
the  scope  this  report.  The  interested  reader  is  referred  to  recent 
literature  addressed  to  this  problem. 

We  have  now  obtained  the  working  equations  for  the  two  approaches  to 
the  problem.  These  are  (7)  and  (10)  for  the  free  energy  minimization  method 
and  (19)  and  (20)  for  the  equilibrium  constant  method.  Both  methods  require 
as  input  the  chemical  potentials  of  the  individual  chemical  species  (see 
equations  (10)  and  (16)).  In  the  ideal  gas  limit  appropriate  to  our  calcula¬ 
tions,  these  and  other  thermodynamic  properties  of  the  individual  species 
can,  in  principle,  be  simply  calculated  by  summation  over  the  species  energy 
levels.  In  practice,  there  can  be  problems,  however,  as  we  shall  mention 
below. 


4.  TRANSFORMATIONS  AMONG  SETS  OF  REFERENCE  SPECIES 

Either  of  the  two  sets  of  equations  obtained,  i.e.  (7)  and  (10)  or  (19) 
and  (20)  can,  in  principle,  be  solved  for  the  concentrations  of  the  various 
constituents  in  the  mixture.  These  equations  contain  the  assumption  that 
there  is  available  a  complete  specification  of  the  system  for  input,  i.e, 
the  temperature  and  density  of  the  mixture,  the  thermodynamic  properties  of 
the  constituents,  and  the  specification  of  the  mixture  Itself  (through  the 
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values  assigned  to  the  x.).  We  have  formally  assumed  that  all  derived  species 
are  obtained  directly  from  the  chemical  combining  of  the  atomic  species. 

Because  of  this,  the  solution  is  specified  in  terms  of  the  concentrations  of 
these  atomic  species.  This  will  become  somewhat  clearer  below.  As  mentioned 
earlier,  the  specification  of  a  precise  set  of  chemical  reactions  from  which 
the  species  are  derived  is  irrelevant  to  the  specification  of  the  thermodynamic 
state  of  the  system.  The  particularization  of  the  reactions  is  needed  only 
for  converting  the  thermodynamic  formalism  into  a  computational  structure.  It 
follows,  therefore,  that  changing  from  one  set  of  reactions  to  another  produces 
no  fundamental  thermodynamic  changes  in  any  of  the  results  obtained. 

Now  a  particular  set  of  chemical  reactions  will  be  convenient,  from  a 
computational  point  of  view,  for  only  a  limited  set  of  thermodynamic  condi¬ 
tions.  Thus,  for  example,  for  chemical  reactions  in  the  system  H20,  H2,  02, 

H,  and  0,  the  use  of  the  atomic  species  H  and  0  as  reference  species,  while 
thermodynamically  correct,  will  not  be  convenient  from  a  computational  point 
of  view  at  temperatures  sufficiently  low  that  little  dissociation  occurs. 

While  equation  (19)  in  terms  of  atomic  reference  species  is  correct  for  the 
reaction 


2H  +  0  Z  H20 

under  all  conditions,  computational  difficulties  occur  at  low  temperatures 
where  the  concentrations  of  the  atomic  species,  ,  become  extremely  small 
and  where  the  equilibrium  constants  (and  hence  K^)  become  extremely  large. 

Despite  the  fact  that  the  obtained  for  H20  at  low  temperatures  by  solving 
(19)  must  be  of  reasonable  magnitude,  it'  is  possible,  In  the  process  of  carry¬ 
ing  out  the  multiplications  required  in  the  right  hand  Bide  of  (19) ,  for 
numbers  to  be  produced  whose  magnitudes  are  outside  the  limits  set  by  the  com¬ 
puter  design.  This  can  often  be  avoided  by  changing  the  order  of  multiplication 
so  that  small  numbers  like  the  C.  alternately  multiply  large  ones  like  or  by 
taking  logarithms.  A  second  problem  Is  associated  with  the  magnitude  orK^. 
While  this  too  Is  a  soluble  problem,  solutions  are  always  artificial  (e.g.,  by 
representing  L  as  a  product  of  factors). 

A  much  more  reasonable  and  more  physical  way  to  solve  both  these  problems 
exists,  namely  the  replacing  of  the  arbitrarily  chosen  set  of  chemical  reac¬ 
tions  used  by  another  set  from  which  reference  species  can  be  defined  whose 
magnitudes  are  within  the  limits  set  by  the  computer.  For  example,  at  low 
temperatures  it  is  obvious  that  H20  should  be  formed  from  molecular  oxygen 
and  molecular  hydrogen.  That  is,  the  computational  framework  should  be  based 
on  H2  and  02  as  reference  species  with  the  H20  reaction  becoming 


H2 


°2  ? 


h2o 


It  should  be  noted  that  the  coefficient  of  the  derived  species  has  been  kept 
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equal  to  unity  to  be  consistent  with  the  definition  of  the  T  matrix. 

In  this  section  we  shall  indicate  how  such  transformations  of  reference 
species  may  be  carried  out.  We  shall  merely  describe  how  the  elements  of  the 
T  matrix  are  to  be  changes,  the  details  required  to  make  the  transformation 
part  of  a  computational  program  being  left  to  the  Interested  reader. 

Equation  (14)  can  be  written 

Hi  "  -I  t..y.  '  i  »  1 . 1  (,22a) 

•j-1  ^  J 


where  j  referes  to  the  reference  species  and  i  to  the  derived  species. 
Equations  (22a)  Include  the  c  identity  equations  for  the  reference  species 
since  these  also  undergo  transformation.  In  matrix  notation  (22a)  can  be 
written 


(22b) 


-y  jy 

with  y *  an  £  element  column  vector  and  y  a  c  element  row  vector.  T  is,  of 
course  an  £  element  row  by  a  c  element  column  matrix. 

In  carrying  out  a  transformation  of  reference  species,  it  is  essential 
that  the  number  of  final  reference  species  be  the  game  as  the  number  of  orig¬ 
inal  reference  species.  In  Brinkley's  terminology^  this  requires  that  the 
dimension  of  the  vector  space  (i.e.  the  rank  of  the  T  matrix)  be  invariant 
under  the  transformation.  This  requires  the  equations  connecting  the  new 
reference  species  with  the  old  to  contain  as  many  equations  as  unknowns,  i.e. 
the  matrix  associated  with  the  transformation  must  be< square.  In  other  words, 
if  double  primes  refer  to  the  new  reference  species,  there  is  a  subset  of 
equations  (22a)  which  can  be  written 


» t 


& 


k.i^ 


» c 


where  the  t,  ,  are  new  the  elements  of  a  square  matrix.  This  can  also  be 
written  in  tile  vector  notation  of  (22b).  Thus 


-*■"  _  -+ 

H  “-Ts  *  U  (22c) 
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where  the  subscript  s  is  meant  to  indicate  the  square  matrix  associated  with 
the  two  set^  of  reference^species.  Since  T  is_|quare,  (22c)  can  formally  be 
solved  for  p  In  terms  of  p"  by  multiplying  Sy  T  from  the  left*  There 
results 


P 


(22d) 


(22b)  can  now  be  writteg.  in  terms  of  the  new  reference  species  by 
substituting  (22d)  for  p  in  (22b).  There  results 

fl 

■+I  1  ■+ 

u ’  =  T  ■  T“ 1  •  u  D  T 

~a  H  -I 


with  the  elements  of  constituting  the  new  coefficient  matrix.  It  should 
be  noted  that  the  elements  associated  with  the  previous  reference  species 
will  no  longer  be  Kronecker  deltas ,  whereas  those  associated  with  the  new 
reference  species  will  become  Kronecker  deltas.  As  will  become  clear  below, 
these  transformations  must  be  accompanied  with  various  changes  in  the 
identitites  of  rows  and  columns  in  the  original  T  matrix. 

The  transformation  is  not  complete  until  the  constants  x.  have  been 
transformed  to  those  appropriate  to  the  new  reference  species.  In  matrix 
form,  (20)  can  be  written 

C  •  T  ='  £ 

Multiplication  by  the  square  array  l”*  from  the  right  yields 

•v  Q 


C  •  T  •  T”1  =  X  •  I”1 

«*  £>  ^  D 

But  T  •  T  ^  =  Tn  so  that  this  becomes 
-  -s  -1 


C  •  I,  -  X  •  i;1 


and  it  is  clear  that  the  transformations  of  the  ^  constitute  the  elements  of 
the  tractor  equation 
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In  order  to  clarify  the  mechanics  of  the  transformation,  let  us  take, 
as  an  example,  the  set  of  reactions  among  the  species  CCL,  CO,  0^,  C,  and  0. 
When  the  atoms  C  and  0  are  taken  as  reference  species,  the  elements  of  the 
mdtrix  are  giv.en  by 

C  0 


Suppose,  now,  that  it  is  desired  to  transform  to  CO  and  as  reference 

species.  The  transformation  matrix  T  and  its  Inverse  are  easily  seen  to  be 
given  by 
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where  the  shift  of  reference  species  is  indicated  by  a  shift  in  the  location 
of  the  diagonal  unit  submatrix  between  the  dashed  lines. 

To  illustrate  the  transformation  on  the  x. >  suppose  there  to  be  1  mole 
of  C02  present  so  that,  for  C  and  0  as  reference  species,  X]_  ■  1>  X2  =  2* 
Thus 

•  (1) 


and,  for  CO  and  02  as  reference  species,  the  values  X2  “  “  are 

appropriate.  2 

Clearly  transformations  of  this  kind  can  be  carried  out  automatically 
within  the  framework  of  any  computer  program  designed  to  solve  the  equations 
for  the  concentrations.  We  shall  not  describe  how  this  might  be  done  in 
practice.  A  word  of  caution  is  in  order,  however,  and  that  is  that,  as  should 
become  clear  in  the  next  section,  since  the  reference  species  are  changed, 
these  transformations  must  always  be  accompanied  with  a  shift  in  the  energy 
differences  (i.e.  heats  of  reaction)  used  to  produce  derived  species  from 
reference  species. 

5.  SOME  PROBLEMS  IN  THE  CALCULATION  OF  IDEAL  GAS  FUNCTIONS  FOR 

THE  INDIVIDUAL  SPECIES 

The  only  microscopic  information  that  has  thus  far  been  included  In  the 
formalism  Is  contained  In  (1)  and  is  merely  the  labeling  of  each  derived 
species  in  terms  of  its  constituent  reference  species.  The  points  of  entry 
for  the  detailed  microscopic  properties  of  individual  constituents  have  al¬ 
ready  been  passed.  These  are  (10)  for  the  free  energy  minimization  method  ■ 
and  (16)  for  the  equilibrium  constant  method.  For  a  pure  component,  the  chem¬ 
ical  potential  is  the  same  as  the  free  energy  per  mole  so  that  In  (9) 


^(T)  =  (G°  -  E°)  + 
RT~  RT  RT 


where  G°  Is  the  ideal  gas  Gibbs  free  energy  for : the  species  and  (E°)  is  the 
reference  energy  at  zero  kelvin.  The  expression  for  the  equilibrium  constant 
(16) ,  becomes 
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lnKi  "  "  ^GT  ~  V  i  +  L  vij  ~  V  j  "  + 

KT  J  RT  RT  J  RT 


Both  methods,  therefore,  require  the  ideal  gas  Gibbs  function  relative  to  the 
energy  of  a  reference  state,  E°  for  the  individual  species. 


Since  only  energy  differences  have  physical  meaning,  the  choice  of  EQ 
values  would  appear  to  be  entirely  arbitrary.  In  a  reacting  system, 
however,  the  difference  between  for  reactant  and  products  is  physically 
meaningful  and  is,  in  fact,  the  reaction  energy  AEQ/RT  ■  (EQ)  i- 

RT  ^  RT 


at  absolute  zero.  For  this  reason,  the  choice  must  be  made  in  a  manner  which 
is  consistent  among  species.  In  most  cases,  taking  E_  to  be  zero  for  the 
reference  species  avoids  .covplic&tl.ons.  Of  course,  all  choices  are  valid 
provided  only  that  they  are  consistent.  It  should  perhaps  be  pointed  out 
again  that  whatever  reference  point  is  taken  for  the  energy  must  be  examined 
carefully  after  any  transformation  to  new  reference  species  is  carried  out. 


The  quantity  (G°  -  E°)/RT  can  be  calculated  for  each  species  using 
statistical  mechanics  and0 the  details  of  the  internal  structure  of  that  species. 
The  Gibbs  free  energy  for  an  ideal  gas  species  is  related  to  the  partition 
function  through  the  relation 


(G 


-  InQ 


where  Q  Is  the  ideal  gas  partition  function.  Q  contains  translational  and 
electronic  contributions  for  all  species  (except  for  the  electrons  for  which 
only  a  translational  part  is  appropriate) .  For  molecular  species  there  are 
additional  contributions  from  vibrational  and  rotational  energies.  A 
discussion  of  methods  for  the  evaluation  of  such  partition  functions  is  out- 
side  the  scope  of  thiB  report  and  the  interested  reader  is  referred  elsewhere 
for  such  details.  In  thiB  report,  we  assume  that  the  Gibbs  free  energies 
have  been  calculated  properly  and  are  available  as  input  to  the  calculation. 

We  shall,  however,  mention  some  problems  associated  with  the  evaluation  of  the 
partition  functions  since  such  problems  Are  not  always  described  in  the 
literature. 

Perhaps  the  most  serious  such  problem  has  to  do  with  the  actual  diver¬ 
gence  of  these  partition  functions.  There  are  an  infinite  number  of  energy 
levels  just  below  the  ionization  limit  of  an  atom.  These  terms,  if  included 
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In  the  summation,  would  cause  the  partition  function  to  diverge.  At  moderate 
temperatures ,which  can  be  as  high  as  several  thousand  kelvins  this  causes  no 
trouble  since  the  contribution  of  successive  energy  levels,  starting  from 
the  lowest,  tends  to  drop  off  rapidly  and  goes  through  a  broad  minimum  starting 
with  a  relatively  small  quantum  number.  The  series  can  therefore  be  treated 
as  an  asymptotic  series  and  cut  almost  anywhere  after  the  contribution  of 
successive  levels  has  become  negligible.  As  the  temperature  is  raised,  however, 
the  breadth  of  this  minimum  narrows  until  ultimately  there  is  no  level  at 
which  the  contributions  of  successive  levels  become  small.  As  soon  as  thi9 
happens,  the  advantage  of  an  obvious  cut-off  point  is  lost  and,  in  fact,  the 
partition  function  begins  to  depend  strongly  on  the  choice  of  cut-off.  This 
is  obviously  a  signal  that  the  theory  on  which  the  calculation  is  based  has 
become  inadequate,  and  this  is  indeed  the  case.  Care  must  be  taken  even  when 
this  problem  does  not  appear  to  exist.  Exact  partition  function  expressions 
are  often  approximated  by  series  which  are  cut-off  in  a  manner  appropriate 
for  low  temperatures.  When  these  series  are  extended  to  high  temperatures, 
they  sometimes  give  convergent,  albeit  wrong,  results  due  to  the  neglect  of  terms 
not  needed  at  low  temperatures  but  needed  at  higher  temperatures .163  These 
neglected  terms  contain  the  divergence  problem. 

An  interesting  way  of  handling  these  divergences  has  been  devised  by 
Woolley*?  and  studied  in  some  detail  by  several  others, 18  in  this  approach, 
sometimes  referred  to  as  physical  cluster  theory,  no  distinction  is  made  be¬ 
tween  a  molecule  whose  constituent  atoms  are  bound  together  by  means  of 
covalent  bonds  and  a  physical  cluster  of  these  atoms  for  separations  and 
conditions  for  which  the  weaker  van  der  Waals  interactions  are  appropriate. 
Reflection  indeed  leads  one  to  the  conclusion  that  any  distinction  which  might 
be  made  between  these  would  indeed  he  artificial. 

The  divergence  which  arises  from  the  summing  of  the  energy  levels  of 
atoms  and  atomic  ions  over  highly  excited  states  of  the  outer  electron  has  been 
considered  in  a  number  of  different  ways  in  the  literature.  Clearly,  the 
representation  of  the  partition  function  of  the  atom  In  the  mixture  as  a  sum 
over  energy  levels  of  the  isolated  atom  is  an  approximation  to  the  actual  many 
body  problem.  At  ordinary  temperatures,  this  approximation  is  reasonable 
since  the  number  of  atoms  in  any  but  the  lowest  lying  energy  levels  is  quite 
small.  Using  semi-classical  arguments,  it  is  easy  to  see  that  the  Bohr  radii 
for  such  levels  are  quite  small  and  that,  therefore,  the  "paths1'  of  the  outer 
electrons  of  different  atoms  do  not  overlap.  As  the  temperature  is  raised, 
the  situation  changes  and  the  Bohr  orbits  become  sufficiently  large  for  there 
to  be  overlap  even  at  ordinary  densities.  The  problem  has  therefore  become  a 
many  body  problem  whose  partition  function  can  no  longer  be  approximated  by 
a  product  of  one  particle  functions . 

19  20  2L 

Several  workers  '  *  nave  devised  density  dependent  cut-off  methods 

which  are  essentially  variations  of  an  approach  due  to  Urey22a,  Fermi22b'and 
also  contained  in  some  unpublished  work  of  Bethe22c.  In  these  approaches,  the 
.energy  levels  are  summed  only  through  those  quantum  numbers  for  which  the  Bohr 
radius  is  less  than  some  function  of  the  average  interparticle  distance.  The 
problem  new  becomes  density  and  temperature  dependent,  the  density  determining 
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the  average  interparticle  distance  and  the  temperature  the  Bohr  radius  of  the 
outer  electrons  as  well  as  the  probability  of  there  being  close  collisions. 

There  are  thus  two  main  approaches  to  the  problem  of  calculating  the 
ideal  gas  partition  functions  of  atoms  and  atomic  ions  depending  on  the 
temperature  range  of  interest  and  the  ionization  energy.  Up  to  moderately  high 
temperatures,  the  sums  can  be  cut-off  at  a  relatively  low  lying  level  based  on 
the  smallness  of  successive  contributions.  At  high  temperatures,  a  cut-off 
can  be  taken  which  depends  on  temperature  and  density  and  which  may  vary  for  a 
given  atomic  species  from  mixture  to  mixture. 

Since  the  quantum  number  takes  on  discrete  values,  it  is  possible  for 
the  summation  to  take  on  discontinuities  as  a  function  of  temperature  and 
density21.  These  discontinuities  occur  where  a  change  occurs  in  the  final 
quantum  number  accepted  under  the  cut-off  criterion.  A  way  to  reduce  this 
problem  has  been  proposed  by  Woolleyl^  (and  independently  by  Gilmore^^) ,  in 
these  approaches,  the  partition  function  is  divided  into  a  sum  of  two  parts. 

One  of  these  is  the  contribution  due  to  states  below  some  quantum  number 
which  is  always  less  than  n  ,  the  density  dependent  cut-off.  The  other  is  due 
to  the  contribution  of  the  £evels  between  n.  and  n  considered  now  to  be  the 
levels  of  one  electron  oh  an  ionized  atom,  xhe  levels  between  n.  and  n  are, 
in  fact,  considered  to  be  hydrogen-like  on  the  ionic  core.  This  methodnas  the 
advantage  of  allowing  one  to  use  tables  of  ideal  ga9  partition  functicms  and  a 
cut-off  criterion  which  depends  on  the  species  only  through  the  nuclear  charge. 

It  is , therefore,  of  considerable  computational  advantage.  Woolley’s  method 
is  designed  particularly  for  treating  the  ionized  gas  as  an  ionic  solution 
within  the  theoretical  framework  set  up  by  Mayer^^a  and  Poirier. 23b 

It  should  be  clear  from  the  preceding  discussion  that  there  are  several 
ways  in  which  one  can  approach  the  divergence  problems  which  arise  in  the  sum 
over  states  associated  with  the  calculation  of  the  ideal  gas  partition  functions 
of  the  individual  species  in  a  chemical  reacting  mixture.  It  is  imperative, 
no  matter  which  of  these  methods  is  used,  that  the  counting  of  states  for  a 
given  species  above  the  dissociation  or  ionization  limit  for  its  constituents 
be  done  in  a  manner  which  is  consistent  with  the  counting  of  the  states  for 
these  constituents  themselves.  If  this  is  not  done  properly,  a  portion  of  the 
phase  space  for  the  mixture  will  have  been  included  at  least  twice  in  the 
statistical  mechanical  development  for  the  mixture  thermodynamic  properties. 

Thus,  one  might  mistakenly  include  the  states  associated  with  the  separated 
atoms  as  a  molecular  state  with  a  van  der  Waals  intermolecular  potential  and 
include  the  same  states  as  free  atoms  with  a  correction  for  nonideality.  A  ^gc 

consistent  way  in  which  this  problem  can  be  avoided  has  been  described  by  Haar. 

6.  THE  THERMODYNAMIC  PROPERTIES  OF  THE  MIXTURE 

Once  the  thermodynamic  functions  for  the  Individual  species  have  been 
calculated  at  the  temperature  of  Interest  and  once  the  energies  AE^/RT  have 
been  chosen,  It  is  possible,  in  principle,  to  solve  for  the  species  concentrations 
using  either  the  method  of,  free  energy  minimization  (minimization  of  (6)  subject 
to  (7))  or  the  equilibrium  constant  method  (solution  of  (19)  subject  to  (18)) 
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and  co  determine  the  concentrations  of  all  species  at  that  temperature  as  a 
function  of  the  most  useful  thermodynamic  variable  (e.g,  pressure,  density, 
etc.).  Given  these  concentrations  and  the  properties  for  the  individual 
species,  one  can  then  calculate  the  thermodynamic  properties  of  the  ideal  gas 
mixture.  In  this  section,  we  shall  indicate  how  this  can 'be  done.  Real  gas 
corrections  to  these  expressions  are  included  below  in  an  appendix. 

According  to  our  definitions  of  concentration  and  compressibility 
factor,  the  compressibility  factor  for  the  ideal  gas  relative  to  that  at 
standard  conditions  is  simply  a  sum  over  species  concentrations,  the  sum 
going  over  both  the  derived  and  reference  species.  Thus, 


Z  - 


P V 
RT 


(23a) 


where,  as  stated  earlier,  C.  is  the  number  of  moles  of  species  i  relative 
to  one  mole  In  a  reference  state  and  Z  is  the  compressibility  factor,  defined 
to  include  the  number  of  moles,  also  with  respect  to  a  reference  9tate,  As 
stated  earlier,  the  choice  of  the  reference  state  is  a  matter  of  convenience 
in  scaling  these  quantities. 

The  internal  energy  can  be  calculated  from 


E 

RT 


l  E 

lA 

1=1 


_i 

RT 


(23b) 


where  Ei  =  (E°  -  and 

where  (E°  -  E°)^  calcuEate^  from  statistical  mechanics. 

The  entropy  is  given  by 

I  =  Z  “  1^1xilnxi  “  ln(  fo>  3  (23c) 

where  P°  is  the  pressure  of  the  thermodynami c  standard  state.  All  other 
thermodynamic  potentials  can  be  written  in  terms  of  combinations  of  (23a, b,c). 
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Among  the  gas  properties  of  interest  are  the  derivatives  of  the 
thermodynamic  potentials.  These  lead  to  such  properties  as  the  specific 
heat  and  sound  velocity.  Such  properties  are  often  obtained  by  numerical 
differentiation  of  tables  of  the  thermodynamic  potentials. 25  xhe  best  of 
such  methods  requires  considerable  caution  in  its  application.  It  is 
necessary  to  be  careful  in  handling  the  end  points  and  care  must  be  taken 
to  ensure  that  the  data  points  are  neither  too  close  together  nor  too  far 
apart.  Numerical  differentiation  is  entirely  unnecessary,  however,  since 
the  derivatives  can  be  obtained  as  the  solutions  of  sets  of  c  linear 
equations  in  c  unknowns ,  c  being  the  number  of  reference  species  in  the 
mixture2c»26, 27#  Several  authors28,29  have  developed  similar  approaches  in 
which,  however,  the  derivatives  are  expressed  in  terms  of  the  solutions  of  £ 
linear  equations  in  £  unknowns ,  2.  being  the  total  number  of  species  in  the 
mixture.  The  difference  between  c  and  £  is  generally  quite  large.  For  high 
temperature  air,  for  example,  l  ~  30  while  c  =  6.  It  should  be  understood, 
however,  that  the  reduction  from  £  equations  to  c  equations  occurs  only  under 
ideal  gas  conditions  (i.e.  y^  =  1) .  When  real  gas  effects  are  included,  it 

is  always  necessary  to  solve  a  set  of  £  equations  in  £  unknowns  for  the  deriva¬ 
tives.  We  shall  proceed  to  derive  the  linear  equations  for  the  derivatives 
within  the  framework  of  the  equilibrium  constant  method.  The  results  obtained 
are,  of  course,  valid  in  either  approach. 

The  derivatives  required  can  in  all  cases  be  written  in  terms  of  the 
derivatives  of  properties  which  are  additive  functions  of  the  properties  of 
the  individual  constituents.  This  means  that  it  is  only  necessary  to  calculate 
the  derivatives  of  functions  of  the  form 


Y 


L“1 


CiYi 


(24) 


til 

where  is  an  arbitrary  ideal  gas  property  for  the  i 


Thus,  for  the  specific  heat  at  constant  volume  one 


constituent . 
has 


(  »  ) 
v  9T  p 


£ 

where  E  =  £  C.E, 

i=l  1 


The  difference  in  the  specific  heats  is  given  by 
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c  -c 

_E _ V 


1  +  1  C  M) 

z  varp 

i  +  £  (— ) 

Z  ^3p; T 


£ 

where  Z  =  £  C . 

i-1 


and  the  ratio  of  specific  heatgis  then  obtainable  from  y 
sound  velocity  can  be  written 


C  /C  .  The 
P  v 


T  ■  ™IZ  [  1  +  2  (f>I  J 
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Useful  in  aerodynamic  calculations  are  the  following 
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The  above  expressions  contain  derivatives  of  additive  functions  as  in 
(24)  which  are  either  with  respect  to  temperature  at  constant  density  or  are 
with  respect  to  density  at  constant  temperature. 

Let  us  first  consider  the  derivatives  at  constant  density,  i.e.,  we 
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3Y 

shall  consider  the  evaluation  of  the  quantity  ' . 


Now,  from  (24) 


(H) 

'‘ST'p 


i) 

^'31  Jp 


y  c  (Ht) 

i  l3T  ;p 

i°l 


(25) 


i 

The  second  term  contains  the  species  ideal  gas  function  (— )  which  can  be 
computed  directly  from  the  partition,  function  of  the  i**1  species.  The  are 
obtained  as  solutions  of  the  equilibrium  equations.  The  second  term  is 
therefore  simply  a  cumulative  sum  involving  known  quantities.  Since  the  Y^ 
which  appear  in  the  first  term  are  also  known,  the  problem  of  evaluating 
(25)  reduces  to  the  determination  of  the  concentration  derivatives.  We  shall 
do  this  for  the  ideal  gas  (i.e.,  y.^  =  1.)  only,  the  extension  to  conditions 
where  this  is  not  appropriate  (i.e.  y.  ?  1.)  being  left  to  the  interested 
reader.  Taking  the  derivative  of  (19  J  yields 
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where  j  runs  over  the  reference  species  only.  There  are  also  the  conservation 
equations  (20) , 
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Since  v 


the  Kronecker  delta  for  reference  species,  this  can  be  written 


l 

X j  “  ^  a 3  ”  1 »  •  ■  • » c  (27) 

J  i-1  1  J 


where,  again,  j  runs  over  the  reference  species  only.  Differentiating  (27) 
and  using  (26) ,  there  results 
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(29) 


is  a  known  quantity,  with  E  the  internal  energy  of  the  ith  species.  The  c 
31nCk 

unknowns,  ( — — —  )  ,  i.e.  the  derivatives  of  the  reference  species,  constitute 

the  unknown  vector  in  (28).  Having  solved  these  equations,  e.g.,  by  inverting 
the  matrix  A,  (26)  can  6e  used  to  calculate  the  derivatives  of  the  ordinary 
species. 

Let  us  now  consider  derivatives  with  respect  to  the  density  at  constant 
temperatures,  i.e,  (3Y)  Differentiation  of  (19)  leads  to  the  relations 
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On  the  other  hand,  differentiation  of  (27)  yields 
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Combining  the  latter  with  (30)  leads  to 
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where 
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(31) 


and  where  A,,  is  again  given  by  (29) .  Since  the  coefficient  matrix  A  ,  is 
identical  in^both  cases,  once  can  use  essentially  the  same  computer  program 
to  solve  (28)  and  (31).  These  solutions  are  then  substituted  into  (26)  and 
(30)  to  obtain  the  derivatives  of  the  remaining  species. 

It  should  be  noted  that  equations  (28)  and  (31)  are  each  c  linear 
equations  in  the  derivatives  of  the  c  reference  species.  At  this  point  it 
should  be  clear  that,  fundamentally,  the  equilibrium  constant  approach  to  the 
problem  is  one  of  solving  c  non-linear  equations  in  c  unknowns. 


As  mentioned  above,  the  ability  to  express  the  derivatives  In  terms  of 
c  equations  in  c  unknowns  is  peculiar  to  the  ide'al  gas  approximation.  As 
soon  as  concentration  dependent  terms  are  introduced  in  the  activity  co¬ 
efficients,  concentrations  of  derived  species  appear  on  the  right  hand  side 
of  (19)  and  hence  on  the  right  hand  sides  of  (26)  and  (30).  This,  then, 
requires  the  solution  of  linear  equation  in  .unknowns  for  the  evaluation 
of  the  concentration  derivatives.  It  may  still  be  possible,  however,  to 
use  these  ideal  gas  expressions  for  the  derivatives  under  conditions  where 
the  concentrations  include  real  gas  effects.  This  can  be  done  when  y ^  4  1. 


but  (  -^-)T  and 


P 


are  negligible  for  major  constituents. 


This  can  even 


be  extended  to  situations  in  which  these  derivatives  are  small  but  not 
negligilbe  by  developing  additive  expressions  for  corrections  due  to  the 
contribution  of  these  derivatives  of  the  y^,  These  expressions  could  then 
be  evaluated  by  an  approximate  numerical  scheme. 


7.  NUMERICAL  METHODS 


A.  INTRODUCTION 

We  now  turn  to  the  central  problem  of  this  report,  namely  that  of 
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solving  for  the  species  concentrations.  This  can  be  done  within  the  frame¬ 
work  of  either  the  free  energy  minimization  method  or  the  equilibrium 
constant  method. 3^  These  are,  of  course,  entirely  equivalent  methods31  and 
the  decision  to  choose  one  over  the  other  is  mainly  a  matter  of  taste.  For 
the  purposes  of  this  report,  numerical  methods  will  mean  those  methods 
appropriate  to  digital  calculators.  It  is  interesting  to  note  that  chemical 
equilibrium  equations  have  been  solved  using  analog  computers.32  Such 
methods  are  considered  to  be  outside  the  province  of  this  report,  however. 


Within  the  framework  of  the  two  methods  the  literature  can  be  further 
divided  between  numerical  approaches  which  are  specifically  tailored  to  a 
particular  chemical  problem  and  those  which  are  general  purpose  approaches. 
Special  purpose  approaches  mainly  predate  the  advent  of  high  speed  electronic 
computers.  At  that  time  the  carrying  out  of  involved  algebraic  operations 
was  much  to  be  preferred  over  long  and  tedious  numerical  calculations  which 
had  to  be  done  on  a  slide  rule  or  desk  calculator.  These  special  purpose 
approaches,  in  the  main,  consist  of  the  reduction  of  the  c  equations  (21)  to 
one  or  two  equations  by  successive  substitution.  With  the  increasing 
availability  of  large  scale  high  speed  computers  whose  main  purpose  is  just 
the  carrying  out  of  tedious  repititious  numerical  operations ,  the  need  for 
the  development  of  such  special  purpose  schemes  has  disappeared. 

Intermediate  between  these  special  purpose  schemes  and  general  purpose 
methods  are  a  number  of  approaches33  which  apply  to  systems  for  which  the 
concentrations  of  those  derived  species  which  contain  more  than  one  reference 
species  is  small.  These  approaches  are  based  on  the  fact  that  the  neglect 
of  such  species  serves  to  uncouple  equations  (20)  thereby  converting  the 
problem  Into  that  of  solving  a  set  of  c  Independent  non-linear  equations. 

Such  approaches,  where  they  are  appropriate,  can  be  useful  in  preliminary 
studies  of  chemical  systems. 


It  Is  interesting  to  note  that,  from  a  numerical  point  of  view,  both 
the  free  energy  minimization  and  the  equilibrium  constant  methods ,  require 
the  determiniation  of  a  set  of  values  of  the  concentrations  of  the  reference 
species  which  minimizes  some  function  of  these  concentrations.  This  is 
obvious  for  the  free  energy  minimization  method  for  which  the  free  energy 
itself  is  being  minimized.  In  the  equilibrium  constant  method,  on  the  other 
hand,  the  exact  solution  of  a  set  of  equations  is  required.  From  a 
numerical  point  of  view,  this  can  also  be  cast  into  a  form  requiring  a 
minimization,  namely  the  minimization  of  an  error  function  which  represents 
departures  from  the  solution. 34  Thus,  the  solution  of  an  arbitrary  set  of 
equations 


i  ■  1,  . . .  ,n 


(32) 


is  equivalent  to  the  minimization  of  the  function 
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(33) 


with  respect  to  the  parameters  x^(i  =  1, 


. ,  n.)  .  This  replacement  by  a 
minimization  problem  is  not  unique,  however.  In  particular,  if  a. .  are  the 
elements  of  any  positive  definite  matrix  a,  then  the  solution  of  the  set  of 
equations  (32)  is  also  equivalent  to  the  minimization  of  the  function 
4>  =  E  C'j  with  both  minima  being  identical.  (33)  is  thus  seen  to  be 


ij 


i  “ij 


the  special  case  where  gt  =  I  the  unit  matrix.  Since  both  methods  lead  to  a 
numerical  minimization  problem,  it  is  possible  to  describe  the  solution  of 
the  equations  appropriate  to  either  in  the  same  numerical  mathematical 
language. 


Non-linear  problems  generally  cannot  be  solved  in  closed  form.  In  the 
present  case,  the  fact  that  the  problem  is  non-linear  can  perhaps  best  be 
9een  in  equations  (21)  of  the  equilibrium  constant  method  where  the  unknowns, 
C,,  appear  as  products  raised  to  powers  (some  of  which  are  negative).  Non¬ 
linear  equations  of  this  kind  will  certainly  not  be  soluble  in  closed  form, 
particularly  in  such  a  way  as  to  be  capable  of  handling  different  chemical 
systems.  It  should  be  noted  in  this  regard,  that  each  chemical  system 
requires  the  use  of  a  different  T  matrix  and  hence  a  different  set  of 
exponents  v  in  (21).  This,  in  turn,  requires  the  solution  of  a  different 
set  of  non-linear  equations.  In  the  absence  of  a  general  closed  form  method 
of  solution,  therefore,  it  is  necessary  to  develop  a  purely  numerical 
approach  to  the  solution  of  the  problem  of  minimizing  <j>  in  (33). 


Such  methods  start  with  a  guess  at  the  set  of  variables  x- . in 

(32).  This  guess  is  substituted  into  4>  in  (33)  and  generally  found  not  to 
be  the  set  of  variables  which  minimizes  <j>.  Procedures  are  then  invoked  to 
improve  on  the  guess  in  a  systematic  way.  If  the  improvement  procedure  is 
properly  designed,  that  set  of  values  of  the  variable  x^,,..^  is  ultimately 
determined  which  does  in  fact  minimize  the  function  if>  within  a  preassigned 
tolerance. 


Clearly,  whether  use  is  made  of  the  free  energy  minimization  method  or 
theequilibrium  constant  method,  and,  in  fact,'  regardless  of  whatever 
numerical  method  is  chosen,  there  must  ultimately  result  a  minimum  value  of  <t>. 
Since  the  function  <j>  describes  a  surface  in  the  n-dimensional  space  whose 
coordinate  axes  are  the  concentrations  of  the  reference  species,  all 
numerical  methods  which  can  be  devised  for  finding  the  minimum  of  <J>  must  be 
geometrically  equivalent  to  starting  from  some  initial  point  on  this  <j>  sur¬ 
face  and  searching  for  the  surface  minimum.  For  this  reason  all  such  methods 
are  referred  to  as  search  methods.  The  function  being  minimized  (i.e.  <(>  in 
(33))  Is  referred  to  as  the  objective  function. 

Non-linear  problems  which  are  fundamentally  difference  from  each  other 
make  use  of  different  objective  functions.  As  will  be  seen  below,  it  can 
also  be  true  that  different  numerical  approaches  to  the  same  non-linear 
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problem  can  be  based  on  the  use  of  different  objective  functions.  In  the 
latter  case,  it  should  be  obvious  that  the  n-dimensional  surfaces  corres¬ 
ponding  to  the  different  objective  functions  should  have  the  same 
corrdinates  (but  not  necessarily  the  same  value)  for  the  objective  function 
at  their  respective  minima.  This  can  usually  be  shown  to  be  trivially  true 
since  these  objective  functions  most  often  differ  from  each  other  by  functions 
which  vanish  at  the  location  of  the  minimum.  Even  when  this  is  not  so,  they, 
of  necessity,  differ  by  functions  whose  surfaces  also  have  minima  at  this 
same  location. 


There  are  a  number  of  objective  functions  which  have  been  widely  used 
in  the  numerical  methods  associated  with  the  direct  free  energy  minimization 
and  equilibrium  constant  methods.  Consider  the  direct  free  energy  minimiz¬ 
ation  method.  Here,  the  objective  function  can  simply  be  taken  to  be  the 
free  energy  of  the  mixture.  In  another  numerical  approach  to  this  method, 
constants  times  the  sum  of  the  squares  of  the  left-hand  sides  of  (7)  are 
added  to  the  free  energy  and  the  total  function  taken  as  the  objective 
function.  Since  the  left-  hand  sides  of  (7)  vanish  at  the  solution,  this 
addition  to  the  free  energy  has  no  effect  on  the  solution.  In  fact,  since 
these  added  quantities  are  positive  away  from  the  solution,  they  tend  to 
increase  the  magnitude  of  the  objective  function  away  from  the  solution  over 
the  value  it  would  have  if  it  were  simply  the  free  energy.  As  a  result,  they 
tend  to  magnify  the  depth  of  the  minimum  in  the  objective  function  over  that 
obtained  fusing  the  free  energy  minimum  alone.  Although  other  objective 
functions  are  possible,  they  will  not  be  included  in  this  discussion. 

While  the  equilibrium  constant  method  can  likewise  be  stated  in  more 
than  one  way  according  to  the  choice  of  objective  function, 36  we  shall  re¬ 
strict  ourselves  to  that  based  directly  on  satisfying  the  relevant  equations. 
Thus,  the  set  of  equations  (21)  will  be  satisfied  if  and  only  if  the  proper 
values  are  used  for  the  reference  species.  In  the  spirit  of  equation  (33) 
it  is  possible  to  define  a  set  of  functions  of  the  concentrations  of  the 
reference  species 


V^(C^, • • • jC^) 


=J1vijKi(a/p<,) 
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The  problem  of  solving  (21)  is  then  equivalent  to  the  minimization  of  the 
objective  function 


(j)  <=>  t 

i"l 


2 


It  should  be  noted  that  the  addition  of  the  Gibbs  free  energy  to  this'  <J> 
results  in  yet  another  possible  objective  function  and,  in  fact,  is  essent¬ 
ially  one  of  the  objective  functions  described  for  the  direct  free  energy 
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minimization  method. 

It  should  be  obvious  that  the  same  search  method  can,  in  principle,  be 
used  for  all  minimization  problems.  For  this  reason  part  of  what  follows 
contains  descriptions  of  and  references  to  general  search  methods  (i.e.  not 
specific  to  the  chemical  equilibrium  problem).  In  practice,  however,  it  is 
often  found  that  search  methods  need  to  be  empirically  tailored  to  the 
particular  non-linear' problem  at  hand.  This  is  not  to  say,  in  the  problem 
under  discussion,  that  the  search  methods  must  necessarily  be  tailored  to  the 
particular  chemical  system  of  interest.  Rather,  these  methods  must  be 
adjusted  to  handle  the  chemical  .equilibrium  problem  in  a  manner  different 
from  an  adjustment  made  for  some  other  non-linear  problem.  This  says,  in 
essence,  that  variations  of  conditions  within  the  same  non-linear  problem 
(e.g.  the  variation  due  to  changing  the  chemical  system  in  the  equilibrium 
problem)  will  result  in  objective  function  surfaces  which  are  much  more 
nearly  alike  than  are  the  surfaces  associated  with  objective  functions  for 
quite  different  non-linear  problems. 

Each  chemical  problem  might  be  expected  to  have  an  optimum  numerical 
method  of  solution.  The  chances  of  there  being  such  an  optimum  numerical 
method  of  solution  for  the  wide  range  of  chemical  problems  associated  with  all 
possible  variations  in  the  T  matrix  is  quite  small,  however.  We  shall,  in 
fact,  assume  that  no  such  general  optimum  method  exists  and  shall  aim  at  the 
description  of  general  search  methods  which  should  have  a  high  probability  of 
converging  to  the  correct  answer  for  almost  all  chemical  systems.  It  should 
be  pointed  out,  however,  that  there  can  also  be  considerable  advantage  to 
making  small  adjustments  on  the  search  method  for  each  particular  chemical 
system  where  such  adjustments  are  feasible.  A  particular  method  with  which 
the  author  has  considerable  experience  and  which  has  been  found  to  converge 
rapidly  in  a  manner  which  is  independent  of  the  initial  guesses  will  be 
discussed  below  in  some  detail. 

A  general  philosophical  point  relating  to  search  methods  needs  to  be 
emphasized.  That  is,  there  are  no  restrictions  on  the  procedures  which  can 
be  devised  for  going  from  point  to  point  on  the  if>  surface,  provided  only 
that  a  proper  test  is  used  for  the  determination  of  when  the  surface  minimum 
is  reached.  In  particular,  the  precise  method  for  the  development  of  the 
path  to  the  minimum  need  have  no  relationship  to  the  chemistry  of  the  problem 
and,  for  that  matter,  none  to  its  mathematics.  In  fact,  it  is  not  even 
necessary  to  follow  a  point  having  a  particular  value  of  $  by  one  with  a 
smaller  value  of  ■,  if  this  turns  out  to  be  convenient.  As  might  be  expected, 
however,  search  methods  have  been  developed  in  which  a  sequency  of  points  is 
traced  out  on  the  surface  in  such  a  way  that  the  $  associated  with  a  given 
point  is  indeed  guaranteed  to  be  less  than  the  s*>  associated  with  the  preceding 
point.  The  ability  to  do  this  guarantees  convergence  to  the  answer  in  the 
sense  that  successively  decreasing  $  ultimately  leads  to  its  minimum  value. 
This  guarantee  is  sometimes  obtained  at  the  cost  of  slow  convergence,  however. 

The  starting  point  has  no  more  a  priori  physical  significance  than  a 
point  used  in  any  other  iteration.  It  follows,  therefore,  that  the 
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particular  method  chosen  for  determining  initial  guesses  need  have  no  relation 
Co  the  chemistry  or  mathematics  of  the  problem  just  as  for  other  iteration 
points.  Many  of  the  methods  as  described  in  the  literature  require  starting 
from  points  which  satisfy  the  mass  balance  equations  (20).  Such  a  restriction 
is  entirely  unnecessary,  however,  since  a  proper  test  for  the  solution  will 
automatically  guarantee  that  the  equations  (20)  are  satisfied  at  the  accepted 
solution.  In  what  follows,  we  shall  have  little  to  say  about  methods  for 
choosing  starting  values. 

While  it  is  generally  not  possible  rigorously  to  prove  convergence  for 
a  search  method  except  in  the  vicinity  of  the  minimum,  empirically  it  is 
found  that  most  search  methods  do  converge  even  from  points  quite  far  removed 
from  the  solution.  The  ease  with  which  this  can  be  done  for  a  particular 
method,  or  equivalently,  the  sensitivity  of  the  method  to  the  choice  of  a 
starting  point,  depends  very  strongly  on  the  method  itself. 

It  should  be  clear  from  the  preceding  discussion  that  central  to  the 
solution  of  non-linear  minimization  problems  is  the  proper  choice  of  a 
criterion  by  means  of  which  one  determines  that  the  solution  has  been  found. 
Criteria  for  solutions  which  are  appropriate  to  the  problems  of  interest 
here  are  obvious.  Despite  this,  they  are  not  always  used.  Clearly,  every 
surface  minimum,  regardless  of  the  problem,  is  characterized  by  the  require¬ 
ment  that  the  derivatives  of  the  objective  function  with  respect  to  the 
unknown  variables  are  each  less  than  some  arbitrarily  small  value.  This 
criterion  would  appear  to  be  a  natural  one  for  the  direct  free  energy 
minimization  method.  In  the  case  of  the  equilibrium  constant  method,  there  is 
the  further  requirement  that  the  equations  (19)  and  (20)  themselves  must  be 
simultaneously  satisfied,  (or  equivalently,  that  <j>  in  (33)  is  sufficiently 
small).  An  intermediate  method  can  also  be  developed  by  using  the  direct 
free  energy  minimization  method  of  solution  but  requiring  that  (13)  be 
satisfied  at  the  solution,  the  being  given  by  (8).3^  Despite  these 
rather  obvious  criteria,  however,  search  methods  all  too  often  use  as  a 
criterion  for  solution  the  condition  that  changes  In  all  the  unknowns  from 
one  iteration  to  the  next  become  vanishingly  small  at  the  minimum. 3^  Although 
the  latter  condition  generally  yields  the  same  results  as  the  former,  it  Is 
not  guaranteed  to  do  so. 

B.  SEARCH  METHODS 

In  recent  years,  an  entire  literature  has  developed  which  deals  with 
non-linear  search  methods. 38  This  literature  consists  mainly  of  descriptions 
of  various  methods  with  occasional  reports  of  experience  with  particular 
problems.  There  is  very  little  in  the  way  of  general  proofs  either  for 
optimization  procedures  or  for  convergence  except,  perhaps,  for  points  near 
the  solution  where  linearization  19  valid.  We  shall  only  describe  the 
general  characteristics  of  non-linear  search  methods  leaving  the  reader  the 
option  of  going  to  the  literature  for  more  detailed  reviews  or  for  the  details 
of  specific  methods. 

The  decision  to  choose  one  particular  search  method  over  all  others  will 
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generally  depend,  at  least  In  part,  on  such  non-mathematical  criteria  as 
relative  speed  of  attaining  the  solution  and  ease  of  programming.  It  will 
also  depend  on  how  often  a  particular  problem  is  to  be  solved  and  for  how 
many  different  chemical  systems.  As  computing  machines  become  faster,  Che 
need  for  a  simply  programmable  method  tends  to  outweigh  considerations  of 
machine  speed.  On  the  other  hand,  the  availability  of  generalized  subprograms 
for  the  complicated  mathematical  and  logical  details  of  a  method  tends  to 
reduce  the  programming  time  required. ^9  it  j.s  clear  from  this,  that  criteria 
for  the  choice  of  one  search  method  over  another  depend  to  a  large  extent,  on 
the  details  of  the  computing  facility  available  for  the  solution  of  the 
problem.  In  the  following,  we  shall  describe  a  number  of  different  search 
methods  in  general  terms  and  shall  include  alternatives  which  might  serve  to 
reduce  computing  time.  The  coverage  of  these  methods  is  not  meant  to  be 
exhaustive.  More  complete  discussions  will  be  found  in  the  several  review 
articles  included  among  the  references  which  deal  with  non-linear  search 
methods . 


The  problem  of  finding  the  minimum  of  a  multidimensional  surface  is  a 
special  case  within  the  general  class  of  problems  of  finding  extreme  values. 
Any  method  developed  for  the  express  purpose  of  finding  minima  of  surfaces 
can  be  applied,  with  at  most  minor  modifications,  to  the  problem  of  finding 
surface  maxima  and  vice  versa.  The  latter  are  more  natural  within  the 
framework  of  the  disciplines  of  economics  and  operations  research.  This 
correspondence  was  specifically  taken  advantage  of  in  an  approach  developed 
at  the  Hand  Corporation  in  which  the  equations  appropriate  to  the  method  of 
the  direct  minimization  of  the  free  energy  were  solved  by  means  of  linear 
programming  techniques  originally  developed  for  the  solution  of  operations 
research  problems.  This  approach  is  discussed  quite  lucidly  In  a  series  of 
Rand  reports  and  will  not  be  considered  here. 


A  non-linear  problem  in  a  set  of  unknowns  can  always  be  converted  to  a 
linear  problem  in  a  set  of  deviations  from  the  current  values  of  the  unknowns. 
In  order  to  do  this  rigorously,  it  is  necessary  to  make  the  assumption  that 
such  deviations  are  negligible  compared  to  the  current  values.  This  can  be 
done  in  0-9)  near  C  th*2  kth  iterate  for  the  reference  species.  By  taking 

'flbMJV  fkV 

-  6Cj  and  neglecting  products  of  the  dC^,  there  results  a 

set  of  linear  equations  for  the  SC.,  These  can  be  solved  for  the  dC  and  a 
new  guess  C  ^+1) obtained.  This  method  is  rigorous  very  near  the  solution 
point.  For^points  far  from  the  solution,  however,  there  is  not  even  the 
guarantee  that  the  corrections  SC,  will  be  In  the  proper  direction.  In 
spite  of  this  the  method  has  been^used  successfully.^ 


Another  method  for  linearizing  non-linear  problems  Is  to  take  a  Taylor 
expansion  about  the  current  values  of  the  unknowns  and  neglect  terms  involving 
products  and  squares  of  deviations.  This  method  differs  from  the  simple  ex¬ 
pansion  In  that  the  first  derivatives  at  the  current  guess  point  appear  in 
the  linear  equations.  This  latter  method  of  linearization  is  essentially 
the  Newton-Raphson  method  for  the  solution  of  sets  of  non-linear  equations. 

We  shall  indicate  several  ways  in  which  this  method  has  been  applied  to  the 
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chemical  equilibrium  problem.  A  common  variation  of  this  approach  makes 
use  of  a  Taylor  expansion  of  the  first  derivatives  about  the  minimum.  In 
this  approach,  the  linear  terms  in  the  deviations  contain  the  second  deri¬ 
vatives  of  the  objective  function. 40  The  latter  approach  is  often  referred 
to  as  the  Gauss-Newton  method. 

Search  methods  for  the  development  of  paths  along  the  surface  to  the 
minimum  can  also  be  referred  to  as  iterative  methods.  This  general  class 
of  numerical  methods,  starts  with  a  guess  value,  which  is  substituted  into 
the  equations  to  be  solved.  On  the  basis  of  the  result  of  this  substitution, 
a  new  value  is  determined  which  then  becomes  the  guess  value  for  a  repeat  of 
this  procedure.  This  process  is  continued  until,  given  a  convergent  method, 
a  solution  is  obtained.  Iterative  methods  can  be  written  in  terms  of  a 
function  of  the  variables  which  gives  a  prescription  for  calculating  the  next 
guess  from  a  given  guess  point.  In  particular,  this  function  defines  a 
sequence  of  operations 


^.<k+l)  +(k+l) 

x.  “  r  Cx  ) 


(34) 


The  functions  T  can  be  viewed  as  operators  which,  when  operating  on  a  guess 
value  x  ,  produce  the  next  guess  value  This  operation  can  also  be  given 

a  geometrical  interpretation  by  writing  the  sequence 


-(Icfl)  -v(k) 
X  ■  X  + 


(35) 


where  p  is  a  unit  vector  in  the  direction  of  the  next  guess  with  X  the  magni¬ 
tude  of  the  step  in  that  direction.-  (35)  is  generally  used  to  describe  descent 
methods ,  pT  then  always  being  taken  as  a  vector  in  the  general  direction  of  the 
minimum. 

Iteration  methods  are  generally  described  as  being  either  of  the  Newton- 
Raphson  type  or  of  the  descent  type,  although  the  distinction  is  not  always 
clear.  In  the  Newton-Raphson  approach,  T(x)  is  obtained  from  a  truncated, 

Taylor  series  about  the  solution.  For  example,  suppose  that  the  surface  minimum 

occurs  at  the  point  x(o^  with  x(1^  a  point  close  by.  It  follows,  then,  that 

^  =  +  J  (x^)  £Co>  -  SWj  where  J(x)  is  the  matrix  of  the 

partial  derivatives  of  v(x),  i.e.  J  "  .  .  If  is  at  the  minimum, 

it  follows  that  \Jj(x(o^)  =  0  so  that,  solving  for  x(°\  there  results 
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*f°*  j-1  <^J  *Cx  P 

If  in  a  sequence  of  Iterations  the  kC^  point  1b  sufficiently  near  th|^ 
minimum,  it  follows  that  the  solution,  can  be  obtained  as  the  (k  +  1)  point 

by  salving  the  equation 


jcw-d.  ;<u_J-i(i<um;<w) 


While  this  is  only  rigorously  correct  at  the  solution,  it  is  used  in  the 
Newton-Raphson  approach  at  all  points  on  the  surface  as  a  means  of  determining 
a  sequence  of  pointB  ,  x(k+2)  ^  aa>|  etc.  which  define  a  path  to  the 

minimum.  It  should  be  noted  that,  for  the  Newton-Raphson  method  the  T 
operator  is  given  by 


-  x-  j"1  <*)  ^(S 


(37) 


Substitution  of  (37)  into  (34)  and  comparison  with  (35)  shows  that  there  is  a 
descent  version  of  the  Newton-Raphson  search  method  with  the  direction  vector 
given  by 


Xv  -  -J-1  CjT)  *(S 


At  the  solution,  (37)  becomes  T(x.  ^  ^  $  ^o  that  x'  x  and  the 

differences  between  successive  iterations  eventually  vanish,  As  has  been 
mentioned,  this  can  be  used  as  a  basis  for  a  criterion  for  solution.  Thus,  a 
solution  is  said  to  have  been  reached  when  the  change  between  successive 
values  of  all  components  of  the  vector  x  become  less  than  some  fraction  of 
their  current  values.  Unfortunately,  this  is  not  a  necessary  and  sufficient 
condition  for  a  solution.  Thus,  while  we  have  shown  this  property  to  hold  at 
the  solution,  (necessary  condition),  it  is  possible  that  for  a  particular 
search  method  a  point  on  the  surface  can  satisfy  this  condition  away  from  the 
solution.  Thus,  satisfying  this  criterion  does  not  guarantee  a  solution 
(sufficient  condition) . 
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Mention  has  been  made  of  a  version  of  the  Newton-Rapheon  method  (often  ^ 
referred  to  as  the  Gauss-Newton  method)  which  makes  use  of  second  derivatives. 
This  can  be  written  in  vector  notation  in  a  form  similar  to  (36).  For 

x;(k+l)  -  x<k)  small,  the  matrix  of  first  derivatives,  J(x)  can  be  expanded 

+lr  +(k) 

about  the  point  x  to  terms  linear  in  the  difference  x  -  x  .  This  results 


-  J(5(k>)  +  H(x(k>)  .  <*<k+1>  -  x(M) 

»X, 

where  H(x)  is  the  matrix  of  second  derivatives  generally  referred  to  as  the 

Hessian  matrix.  If  the  assumption  is  made  that  x  i9  the  solution  point, 

-*(k+l) 

it  follows  that  J(x  )  =  0  and  the  above  relation  becomes 


0  «  J(x(k))  +  H(x^k)) 


(;<k+1>  -  j<k)) 


or: 


H(J(W)  -  H(I(k+1))J(k>  -  J<S(k)) 


from  which 


-  H'1  (J<k)) 


quite  analogous  to  (36).  In  this  case  (x)  is  given  by 


T(x)  =  x  -  H  1(x)  J  (x) 

In  what  follows  we  shall  discuss  the  Newton-Raphson  method  only  In  the 
form  (36) . 


Many  variations  of  the  Newton-Raphson  search  algorithm  (i.e.  (36)  have 
been  devised  through  variation  of  the  definition  of  J(x).  These  alternate 
methods  are  mainly  used  when  the  calculation  of  J (x)  and  its  inverse  either 
involves  excessive  computing  time  or  is  overly  complex.  The  simplest  varia¬ 
tion  involves  a  simple  iteration.  In  that  case  J(£)  =  XI,  with  I  the  unit  matrix, 
and  the  algorithm  for  the  choice  of  successive  values 

T(x)  -  x  -  (x)  AiKx) 

which,  at  the  solution  becomes  T(x,  .)  =  x,  .  As  will  be  seen  below  this 
approach  falls  within  a  category  oF  methods  which  we  have  designated  direct 
search  methods.  An  appropriate  value  of  A  can  be  found  empirically. generally 
through  monitoring  the  behavior  of  the  x  vector  in  successive  Iterations.  The 
value  of  A  used  will  depend  strongly  on  the  nature  of  the  function,  ^(x) .  For 
example,  consider  the  one-dimensional  problem  ty(x)  =  0.  If  ^(x)  is  such  that 

<  x(o>  implies  $(x)  >  0  then  clearly  A  >  0  is  required.  If,  -on  the 

other  hand,  x^>  x^°^  implies  i|j(x)  <  0,  A  <  0  is  required.  These  require¬ 
ments  would  easily  be  obtained  empirically. 
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A  more  complicated  variation  of  the  Newton-Raphson  method  involves  re¬ 
placing  the  "tangent"  matrix  J(x)  by  a  "secant"  matrix  In  this 

approach,  the  partial  derivatives  in  the  J  matrix  are  replaced  by  finite 
differences  based  on  neighboring  solution  points.  For  an  n  dimensional  $ 
space  (i.e.  for  n  unknowns),  the  secant  approximation  requires  having  in  hand 
n  +  1  neighboring  values.  ^The  method  thus  does  not  really  get  under  way 
until  the  objective  function  has  been  evaluated  at  n  +  1  points.  The 
procedure  can  be  started  by  choosing,  somewhat  arbitrarily,  these  n  +  1 
points,  evaluating  the  function,  calculating  the  elements  of  the  matrix  J 
and  proceeding  to  invert  J  as  in  the  Newton-Raphson  method.  Since  the 
derivatives  of  the  function  do  not  have  to  be  evaluated,  computing  time  is 
reduced.  Since  the  secants  make  use  of  stored  evaluations  of  the  objective 
function  at  the  n  +  1  points,  the  programming  is  also  simplified.  Since  the 
method  still  requires  a  matrix  inversion  for  determining  J  (x)  from  J(x), 
the  saving  in  machine  time  will  be  modest.  A  dramatic  reduction  in  computational 
time  can  be  had  by  means  of  a  variation  on  the  secant  version  of  the  Newton- 
Raphson  method  in  which  each  component  x  of  the  unknown  vector  ^  is  treated 
separately  thus  eliminating  the  need  for^matrix  inversion.  That  approach  is 
associated  with  the  name  of  Wegstein^  antj  will  be  described  below  among  the 
direct  search  methods. 


An  interesting  variation  of  the  Newton-Raphson  method  has  recently  been 
described. 41  In  that  approach,  the  coordinate  system  in  parameter  space  is 
rotated  into  the  space  spanned  by  the  eigenvectors  of  J,  In  that  space,  J 
can  be  written  in  diagonal  form.  The  advantage  of  this  approach  is  that  since 
the  new  "parameters"  (i.e.,  the  eigenvectors  of  J)  are  orthogonal,  they  do  not 
interact  with  each  other  so  that  the  minimization  can  be  carried  out 
lndpendently  for  each. 


A  number  of  methods  have  been  developed  for  the  solution  of  chemical 
equilibrium  problems  which  make  use  of  a  Newton-Raphson  search  method.  The  ^ 
two  most  widely  used  are  one  due  to  Brinkley^  and  a  second  due  to  Huff  et  al. 
The  relationship  ,between  them  is  discussed  in  some  detail  In  the  review  of 
Zelesnik  and  Gordon.^  Each  method  is  sufficiently  well  described  In  the 
literature  to  preclude  the  need  for  a  description  here.  Brinkley's  approach 
is,  of  course,  built  around  the  solution  of  a  set  of  equations  whose  number 
is  the  same  as  the  number  of  reference  species.  It  is  interesting  to  note 
that  the  method  of  Huff  et  al,  as  originally  described,  considers  both  the 
reference  species  and  the  derived  species  on  an  equal  footing  thus  requiring 
the  solution  of  £  equation  for  the  £  composition  variables.  Since  £  can 
become  prohibitively  large  for  matrix  techniques ,  methods  for  reducing  this 
requirement  are  of  interest.  Zelesnik  and  Gordon  indicate  that,  as  expected, 
the  method  can  be  converted  to  one  in  which  only  the  reference  species  are 
required,  thereby  reducing  the  problem  to  that  of  solving  c  equations  in  c 
unknowns,  where  c  <  £  is  always  the  case.  They  Indicate,  however,  that  this 
can  result  in  convergence  problems  unless  the  reference  species  are  carefully 
chosen  so  as  to  be  major  constituents  or,  if  not,  unless  the  initial  guesses 
are  good.  This  problem  does  not  appear  to  arise  when  all  £  equations  are 
solved  presumably  since  these  always  contain  the  major  constituents. 
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C.  DIRECT  SEARCH  METHODS 

Search  methods  fall  into  two  categories  on  the  basis  of  the  ease  with 
which  they  can  be  applied  to  new  problems.  This  can  be  of  considerable 
practical  importance  in  the  chemical  equilibrium  problem.  For  example,  after 
programming  an  ideal  gas  calculation,  one  might  wish  to  add  in  real  gas 
effects  by  relaxing  the  restriction  y 1.  For  some  search  methods  this  can 
pose  a  major  reprogramming  problem,  wftile  for  others.,  th.e  conversion  is 
rather  simple.  These  two  classes  of  search  methods  can  be  described,  roughly, 
in  terms  of  the  extent  to  which  they  depend  on  the  content  of  the  equilibrium 
equations . 

Search  methods  of  the  general  Newton-Raphson  type  require  considerable 
detail  since  they  call  for  the  calculation  of  first  (and  sometimes  second) 
derivatives  of  the  objective  function  in  addition  to  the  evaluation  of  the 
objective  function  itself.  Reprogramming  for  such  methods  can  call  for 
considerable  effort.  On  the  other  hand,  search  methods  exist  which  at  each 
step  require  only  that  the  objective  function  be  evaluated.  Such  methods  we 
shall  designate  direct  search  methods  since,  in  the  case  of  the  most  obvious 
one,  the  minimum  in  the  objective  function  surface  is  sought  be  a  direct 
search  in  variable  space.  The  usual  definition  of  direct  search  methods  is 
much  more  narrow  than  this.  We  shall,  in  fact,  include  among  our  direct 
search  methods  certain  ones  that  are  variations  of  approaches  not  normally 
defined  as  being  of  the  direct  search  type,  e.g.  variations  of  the  Newton- 
Raphson  method.  This  underlines  the  arbitrariness  of  these  classifications. 

There  are  several  advantages  to  be  had  from  the  use  of  direct  search 
methods.  Since  only  the  objective  function  is  evaluated,  a  minimum  of  re¬ 
programming  effort  is  required  when  the  objective  function  Is  altered. 
Furthermore,  the  simplicity  of  the  search  aspect  allows  for  simple  initial 
programming  thereby  reducing  the  time  required  for  going  from  the  develop¬ 
ment  of  the  problem  to  the  working  computer  program.  Experience  has  also 
shown  that  direct  search  methods  can  be  successful  for  problems  which  are 
poorly  conditioned  for  the  other  methods .43,44  ^  exaTnpie  might  be  a  problem 

for  which  the  matrix  of  the  derivatives  of  the  objective  function  with 
respect  to  the  unknowns  is  ill-conditioned  in  some  part  of  variable  space 
even  far  removed  from  the  solution.  Since  that  matrix  plays  a  key  role  In 
the  Newton-Raphson  method,  there  can  be  trouble  with  it  If  an  interaction 
happens  to  come  near  such  regions  of  ill-conditioning. 

Perhaps  the  simplest  of  the  direct  search  methods  is  one  that  merely 
involves  testing  the  objective  function  in  a  stepwise  variation  of  the  un¬ 
known  parameters  (perhaps  on  a  grid} ,  accepting  only  those  changes  which 
reduce  that  function.  The  pattern  search  method  of  Hooke  and  Jeeves^  i8 
an  improvement  over  this  simplistic  approach.  Their  approach  makes  use  of 
information  obtained  but  usually  ignored  in  a  simple  step  wise  search  on  a 
grid.  In  their  direct  search  method,  two  kinds  of  multi-dimensional  steps 
on  the  $  surface  are  defined-an  exploratory  step  which,  when  successful 
is  followed  by  a  pattern  step.  For  the  former,  the  unknown  parameters  are 
varied  individually  while  for  the  latter,  the  parameters  are  all  changed 
simultaneously.  Exploratory  steps  proceed  as  follows.  The  objective  function 
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.(o) 


the  parameter  x£°^is  increased  to 


v2  *  * ' 
,  i  ■  e  ■  A.^ , 


(o)v 
X  ) 

n 


is  evaluated  at  the  location  of  the  initial  guess  ,  x, 

and  the  result  stored.  Using  a  prescribed  step  size  for 

+  A^  and  the  objective  function  evalu¬ 
ated.  Should  this  lead  to  a  reduction  in  the  objective  function,  x5°^is 


replaced  by  x’  =»  x^°^+  A.  80  that  the  new  guess  point  becomes  (x! , 
and  the  9ame  procedure  is  followed  for  x~  ,  using  A«,  the  prescribed  step 
size  for  x„.  If,  on  the  other  hand,  replacing  x}  by  x.  +  A  does  not  reduce 

2  *<°>  -  Vs 


x(o)  x(o)) 
2  n  ; 


the  objective  function,  a  change  of  direction  is  attempted,  i.e.  ^ 
tried  in  place  of  x^°^  When  this  too  does  not  lead  to  a  reduction 


in  <J>,  A^  is 


replaced  by  aA- (a<l)  and  the  process  is  repeated  from  the  start.  This  is 
continued,  eacn  new  A_  being  decreased  when  there  is  no  success,  until  there 
is  either  a  decrease  In  the  objective  function  for  some  value  of  A^  or  until  A^ 
becomes  smaller  than  a  preassigned  limit.  In  either  case,  the  same  process  is 
carried  over  to  x„.  In  the  former  case  x  is  changed  by  the  current  value  of 
in  the  direction  indicated  by  the  decrease  in  $  while  in  the  latter  case  x£0' 

is  used  as  x'  and  the  guess  point  is  therefore  not  changed.  This  process  is 
applied  to  each  parameter  in  turn  until  each  of  x^,...,  x  hais  been  varied. 

There  are  now  two  possibilities  -  either  the  new  guess  point  x-  -  (x'  ...,  x*) 

:(o)_ 


=  (x^-,  ...,x^)  or  these  two  points 
In  the  former  case,  the  surface  minimum  has 


is  identical  with  the  initial  guess  point  x 

differ  in  at  least  one  parameter 

probably  been  reached.  At  this  point  a  proper  test  for  solution  is  made  to 

see  if,  indeed,  the  surface  minimum  has  been  reached.  When  x-  differs  from 

-+(o)  ->(1)  -*-(o)  _  -> 

x  in  at  least  one  parameter  a  change  vector  x  -  x 


=  A  has  been  deter- 


mined  which,  when  drawn  from  the  initial  guess  point  xk  points  in  the  general 
direction  of  the  minimum.  Hooke  and  Jeeves  suggest  that,  on  the  average,  this 
change  vector  indicates  the  direction  to  the  minimum  sufficiently  well  enough 
for  the  objective  function  to  be  still  smaller  at  5P  ^  ^  '+  t0-  The  step 

x^=  x(1)+  A  is  called  a  pattern  step.  Of  particular  importance  is  the  fact 
that  all  parameters  are  varied  simultaneously  in  a  pattern  step  whereas  in  an 
exploratory  step  each  of. the  n  parameters  is  varied  at  least  once.  The  time 
required  for  a  complete  exploratory  step  for  n  unknowns  can  obviously  be  much 
more  than  n  times  than  required  for  a  pattern  step.  Hooke  and  Jeeves  choose 
to  restrict  themselves  to  one  pattern  step  after  a  successful  exploratory 
step.  They  then  follow  this  single  pattern  step  by  a  new  exploratory  step. 
Clearly,  one  might  modify  their  method  by  making  the  number  of  consecutive 
pattern  steps  after  a  successful  exploratory  step  a  variable  or  even  by  continu¬ 
ing  with  pattern  steps  until  there  ceases  to  be  a  reduction  in  the  objective 
function. 


A  number  of  further  variants  of  the  pattern  search  method  are  possible. 

One  might,  for  example,  choose  to  vary  each  component  step  of  the  exploratory 

phase  from  the  original  guess  point  x  ^so  that  x^  =  x^+  A  where  A  = 

(Ai*  A^,  ...  A^)  as  calculated.  This  does  not  hold  for  the  method  as  we  have 

outlined  it  above  since  each  component  is  changed  on  a  successful  decrease  in 
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Thus,  in  other  words. 


A  7*  A. 
P 


One  might,  furthermore,  devise  variations  on 


the  pattern  step  itself  by  replacing  the  vector  A  as  determined  by  the  explora¬ 
tory  step  by  a  different  change  vector  t  for  usepin  the  pattern  step.  One 
might,  for  example  determine  as  some  £ind  of  average  over  the  vectors 

used  in  several  of  the  preceding  pattern  steps  or  one  might  simply  take  a 
fraction  of  J  . 


Another  approach  which,  according  to  our  definition,  can  be  included 
among  the  direct  search  methods  is  actually  an  independent  parameter  version 
of  a  variation  on  the  secant  approximation  to  the  Kewton-Raphson  method.  The 
method  is  also  referred  to  as  the  Wegstein  method.  The  basic  algorithm 

for  one  dimension,  is  available  among  the  basic  library  routines  at  many 
computer  installations.  The  method  has  been  used  extensively  by  the  author 
and  his  collaborators  in  the  production  of  thermodynamic  tables  of  atmospheric 
gases. This  method  has  been  found  to  produce  rapid  convergence  even  when  the_^ 
initial  guess  values  differ  from  the  answers  by  more  than  factors  of  10+^  or  10 
When  coupled  with  a  method  for  automatically  converting  from  one  set  of 
reference  species  to  another,  where  necessary,  the  method  should  produce  rapid 
convergence  entirely  independent  of  initial  guesses  and  should  therefore  be 
a  truly  automatic  approach.  We  shall  describe  the  method  in  some  detail  in 
one  dimension  but  shall  not  specifically  relate  it  to  the  actual  Newton-Raphson 
method.  The  extension  to  multi-dimensional  problems  in  the  independent 
parameter  approximation  will  then  be  shown  tD  be  straightforward. 

This  method  can  best  be  described  in  terms  of  the  iteration  (34) .  It 
should  be  noted  that  (34)  can  also  be  applied  to  the  solution  of  equations 
of  the  form  iji(x)  =  0  since  one  can  take  f(x)  in  the  form 


«;«) 


;«♦ 


where  a  is  non-zero.  With  this  form  for  T(x),  (34)  reduces  to  the  identity 
-*(k+l)  -+(k)  ->(k)  _ 

x  =*  x  at  the  solution  where  ^(x  )  ■  0.  In  the  Wegstein  method,  which 

we  shall  describe  in  one  dimension,  (34)  is  modified  by  superimposing  on  it  an 
in-aut  averaging  algorithm.  Thus,  rather  than  take  x^+  ^  as  given  by  (34) 
for  the  next  guess,  an  average  between  the  value  x^into  the  iteration  and 
x^  +  ^  out  of  the  iteration  is  taken  to  obtain  as  the  next  guess 


“(k+1)  .  x(k+l)  _  q(x(H-l)_  xOO)  (38) 

This  produces  a  value  given  by  x^+l)  ag  caiculated  in  (34)  minus  a  fraction  q 

(k)  (k+1) 

of  the  calculated  change  between  x'  y  and  x  .  It  is  easy  to  show  that,  in 

the  secant  approximation,  q  =  ---_  where  m  is  the  slope  of  the  secant  between 
(k-1)  (k)  ra~i 

x  and  x  .  This  enables  q  to  be  calculated  at  each  iteration  step.  We 

are,  however,  much  more  interested  in  the  advantages  of  using  q  as  a  fixed 
empirical  quantity. 

The  use  of  (34)  by  itself  will  lead  to  one  of  the  following  behavior 
patterns  for  successive  values  of  x(^),  (i.e.  as  functions  of  k  the  iteration 
number) , 
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(k) 

1.  The  values  of  x  oscillate  and  converge. 

.  (k) 

2.  The  values  of  x  oscillate  and  diverge. 

(k) 

3.  The  values  of  x  converge  monotonically. 

(k) 

4.  The  values  of  x  diverge  monotonically. 


It  is  also  possible  for  the  character  of  the  behavior  pattern  to  change  as  n 
progresses.  Generally,  however,  one  particular  kind  of  behavior  should 
predominate.  Each  of  the  above  behavior  patterns  for  (34)  can  be  improved 
through  the  use  of  (38)  with  fixed  q  in  the  sense  that  convergence  can  be 
speeded  up.  This  can  be  done,  in  each  case,  by  means  of  q  values  in  the 
following  ranges. 


1.  Oscillatory  convergence 

2.  Oscillatory  divergence 

3.  Monotonic  convergence 

4.  Mono tonic  divergence 


0  <  q  <  0. 5 
.5  <  q  <  1. 
q  <  0 
q  >  1 


The  reason  ,  in  each  case,  follows  from  the  use  of  q  in  (38)  as  an  in-out 
averaging  parameter.  Thus,  for  the  case  where  successive  values  oscillate 
Wlth  toward  the  answer  superimposed  on  this  oscillation,  it  is  clear 

that  x  +~  will  11a  between  x  *  and  x^  .  Thig  is  what  is  meant  by . . 
oscillation.  Each  value  will,  however,  on  the  average  be  closer  to  k  * 

than  to  x  .  This,  in  turn,  is  what  is  meant  bv  convergence.  Thus  q  must 
be  so  chosen  that  x^+l)iies  between  x'k'->  and  x^k+l)  closer  to  x(k+l)  , 
than  to  x<k>.  This,  in  turn,  is  what  is  meant  by  convergence.  Thus  q  must 
ao  chosen  that  x^+l)  lies  between  x^  and  but  closer  to  x^k+^\ 

This  is  clearly  the  case  only  for  0  <  q  <  0.5.  Examination  of  the  other 

ranges  listed  for  q  values,  shows  them  also  to  be  designed  to  produce  values 
of  x(k+l)  in  the  general  direction  of  the  solution  from  Wegstein, 

with  several  examples,  illustrates  the  damping  effect  this  method  has  on 
oscillations,  its  ability  to  change  divergent  behavior  into  convergent  behavior, 
and  even  to  speed  up  convergence  where  there  is  mono tonic  convergence.  It 
should  he  noted  that  q  =  1.0  causes  the  iteration  averaging  scheme  (38)  to 
become  x(k*-l)=  x^k^  which,  destroys  the  ability  to  change  x  from  one  Iteration 
to  the  next  and  hence  causes  a  false  lacking  in  on  an  answer. 


The  extension  of  the  Wegstein  method  to  multi -dimensional  problems  is 
particularly  simple.  It  is  merely  necessary  to  state  such  problems  in  such 
a  way  that  the  search  method  can  be  independently  applied  to  each  variable. 
In  an  approach  used  extensively  by  the  author  and  his  collaborators^,  n 
equations  (38)  were  written,  one  being  written  for  each  variable  with  each 
equation  having  Its  own  q  parameter.'  In  this  approach,  the  algorithm 
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was  taken  for  (34)  for  all  species  except  the  electrons,  the  superscripts 
referring  to  iteration  number.  A  different  scheme  had  to  be  used  for  the 
electrons  since  for  them  Xj  “  0.  The  algorithm 

c  -  -  Y  ,  C(k)  i  -  c  (39b) 
1  i-1  13  1 


was  taken  for  the  electrons,  the  sum  extending  over  all  species  except  the 
electrons  themselves .  It  should  be  noted  that,  according  to  our  definitions, 
t  for  j*c  is  positive  for  negative  ions  and  negative  for  positive  ions  so 
that  the  former  substract  from  the  sum  while  the  latter  add  to  it  as  required. 
Also  note  that  (39b)  is  merely  (20)  for  the  electrons  with  Xj=0» 

The  iteration  process  proceeds  as  follows.  Initial  guesses  for  the 
concentrations  of  the  reference  species  are  substituted  into  (19)  thereby 
producing  initial  guesses  for  the  derived  species.  The  method  has  been 
found  to  be  sufficiently  independent  of  these  initial  guesses  for  the 
reference  species  to  allow  the  same  set  of  initial  guesses  to  be  used  for 
all  problems.  The  initial  guesses  thus  calculated  for  the  derived  species 
are  substituted  into  (39a)  and  (39b)  yielding  interim  new  guesses  for  the 
reference  species.  These,  are  then  substituted  into  (38)  and  new  guesses 
obtained  for  the  reference  species.  These  are  used  in  (19)  to  start  the  next 
Iteration.  This  process  is  continued  until  the  proper  criterion  for  solution 
is  satisfied.  It  should  be  noted  that,  according  to  (20) 


l 

l 

i=l 


1 


at  the  solution  in  which  case  (39a)  becomes  C 


(k+1)  _  (k)  ,  . 

=  C,  as  required. 
J  i 


1 


For  the  various  systems  studied  by  the  author  and  his  collaborators, 
the  search  method  was  found  to  be  insensitive  to  the  choice  of  each  of  the 
q.  except  for  a  small  sensitivity  associated  with  the  q  of  the  electron 
concentration.  This  lack  of  sensitivity,  particularly  for  the  atomic 
species,  persisted  over  an  extremely  large  range  of  conditions.  These 
included  temperatures  and  densities  for  which  the  reference  species  were 
present  in  trace  amounts  with  the  major  constituents  being  molecular  species, 
those  for  which  the  reference  species  dominated,  as  well  as  those  for  which 
the  reference  species  were  again  present' in  trace  amounts  with  atomic  ions 
being  the  main  species.  In  all  such  cases  the  value  q  =0.5  was  used  for 
all  atomic  reference  species.  ^ 


The  electron  concentration  behaved  somewhat  differently  as  a  function  of 
iteration  number,  as  might  have  been  expected  since  (39b)  is  of  quite  a 
different  form  than  is  (39a).  For  the  electron  concentration,  there  tended 
to  be  oscillatory  divergence  with  q  3  0.8  being  required.  Somewhat  more 
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rapid  convergence  was  obtained  by  starting  each  problem  with  q  ■*  0.5  for  the 
electron  concentration  and  increasing  the  value  of  q  slightly  with  each 
iteration  so  as  to  reach  q  *  0.8  after  a  small  nuttber  of  iterations.  It  is 
important  to  point  out  that  this  scheme  (i.e.,  q  =  0.5  for  all  specieB  except 
the  electrons  and  a  slowly  increasing  q  value  for  the  latter)  once  adopted 
was  never  modified  regardless  of  the  chemical  system  studied. 

A  number  of  precautions,  common  to  Iteration  methods  had  to  be  taken. 

A  ceiling  was  placed  on  the  relative  magnitude  by  which  the  C.  were  allowed 
to  change  in  one  iteration.  Thus ,  for  J 


|c.(t+1>  -  c.<k>| 

^ -  >  100  or  <  the  algorithms 


,oo 

'  i 


c  <k> 

C  Ck+D  .  100  c  <k>  and  c  Ocfl)  „ 

j  i  j  loo 


were  used. 


In  practice  these  restrictions  were  found  to  operate  only  for  current  guess 
points  quite  far  removed  from  the  solution.  These  ceilings  in  effect  guide 
the  iteration  point  into  the  neighborhood  of  the  solution.  They  play  an 
important  role  in  making  the  problem  independent  of  the  initial  guess  values 
as  they  undoubtedly  aslo  would  if  applied  to  other  search  methods. 
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APPENDIX- I 

ILLUSTRATIONS  OF  THE  MATRICES  ASSOCIATED  WITH  THE  STOICHIOMETRY  OF  REACTING 
GASEOUS  MIXTURES 


A  number  of  matrices  having  to  do  with  the  stoichiometry  of  coupled 
chemical  reactions  in  gaseous  mixtures  have  been  defined  in  the  text.  A 
complet!  understanding  of  these  matrices  and  of  their  relationships  to 
the  chemistry  of  coupled  chemical  reactions  is  absolutely  essential  here. 
Without  such  an  understanding,  much  of  what  appears  in  thi9  chapter  becomes 
unintelligible.  For  this  reason,  we  include  in  this  Appendix  the  matrices 
associated  with  a  typical  set  of  coupled  reactions.  We  have  purposely  chosen 
a  set  of  reactions  which  include  electron  attachment  and  detachment  (i.e., 
ionization)  in  order  to  illustrate  our  treatment  of  the  electron  as  a 
reference  species.  This  treatment  differs,  in  some  respects,  from  standard 
chemical  treatments,  particularly  with  regard  to  notation. 


The  set  of  reactions  with  which  we  shall  be  concerned  in  this 
Appendix  are  reactions  among  the  following  thirteen  species:  N^,  0^,  ^0, 

NO^,  N+,  0+,  02+,  ’  N++*  0++»  °>  and  e  *  where  e  refers  to  the 

electron,  Equation  &L)  of  the  text  requires  these  to  be  vjritten  N2,  C>2, 

N20,  N02,  Ne_“  ,  0e~^,  °2eI1>  °2e+1»  Nel2’  0e-2’  N>  °»  and  e  >  where  a 


negative  subscript  for  the  electron  indicates  the  absence  of  an  electron 
and  a  positive  subscript  the  presence  of  an  extra  electron.  It  should  be 
noted  that  our  notation  differs  from  standard  notation  only  in  that  the 
usual  practice  of  writing  ionic  charge  as  a  superscript  on  the  chemical 
symbol  has  been  replaced  by  specific  reference  to  the  presence  or  absence 
of  the  electron  treated  as  a  chemical  species. 


This  notation  leads  to  the  following  v 


ij 


matrix  elements. 


i 
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REFERENCE  SPECIES 


N 

0 

e 

Species 

N2 

2 

°2 

2 

n2o 

2 

1 

no2 

1 

2 

Ne”x  (N+) 

1 

-1 

Oe^  (0+) 

I 

-I 

•f  <NJ 

o 

w 

H 

1  1 

Q) 

CN 

o 

2 

-1 

°2<1  (0? 

2 

1 

Ne"2  (N4^) 

I 

-2 

0e“2  <0++) 

1 

-.2 

N 

1 

0 

1 
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It  should  be  pointed  out  that  this  set  of  v . .  matrix  elements  is 
associated  with  the  chemical  symbols  of  the  varioui  species  and  so  is 
permanently  associated  with  these  species,  regardless  of  the  definition  of 
reference  species.  The  elements  of  the  various  T  matrices,  on  the  other 
hand,  depend  very  strongly  on  the  choice  of  reference  species. 

According  to  usual  practice  and  notation,  the  chemical  reactions  to 
which  we  shall  restrict  this  discussion  might' be  written,  with  atoms  and 
the  electron  as  reference  species 


2N  t  N2 

20  2  °2 
2N  +  0  +  N20 

N  +  20  +  K02 

N  t  N++e“ 

■+■  +  — 

0  +  0  +e 

20  +  02+e" 

20  +  e  0  2~ 

N  ?  N+++2e*' 

0  T  0+++2e” 

N  i  N 

0*0 

—  +  “ 
e  ■*-  e 


The  T  matrix  associated  with  equation  (5)  of  the  text  was  a  square 
matrix  and  included  both  a  row  and  a  column  for  each  species.  The  reactions 
above  define  such  a  matrix  as  follows. 
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K 


K 


0, 


V 

no2 

N+ 


++ 


Og  n20  no2  n+  0+  Og+ 


N 


++ 


1 


N  0  e 

-2 

-2 

-2  -1 
"1  -2 
-1  1 

“1  1 
-2  1 
"2  "I 
-1  -2 

"1  -2 

1 

1 

1 


where  a  blank  position  Is  used  to  indicate  a  aero  matrix  elements. 


Somewhat  better  symmetry  can  be  obtained,  in  the  sense  that  reference 
species  always  appear  on  the  same  side  of  the  equation,  if  the  ionization 
reactions  are  written 

N  -  e”  t  N+ 


.0  -  e  +0 


°2  “  e  *  °2 


N  -  2e~+  N4"* 
0  -  2e-*  0++ 


It  should  be  noted  that  this  has  no  effect  on  the  matrix  element  t... 
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■The  second  form  of  the  T  matrix  defined  in  the  text  uses  the  fact  that 
the  column  under  a  derived  species  contains  an  entry  only  for  the  row  assoc¬ 
iated  with  that  same  derived  species,  with  that  entry  always  being  +1.  With 
this  +1  understood  (and,  more  important,  supplied  when  needed)  the  T  maxtrix 
can  be  written  in  terms  of  the  reference  species  only.  This  results  in  the 
following  matrix. 


X 

0 

e~ 

»2 

-2 

_  2 

n2o 

-2 

-1 

no2 

-1 

-2 

u+ 

-1 

1 

0+ 

-1 

1 

°2+ 

-2 

1 

-2 

-X 

irH‘ 

-1 

2 

o*+ 

-1 

2 

N 

1 

0 

1 

e~  1 


It  should  be  noted  that,  because  atoms  and  the  electron  are  used  as  reference 
species,  these  matrix  elements  correspond  exactly  to  the  v  .  That  this  is 
not  so  for  other  choices  of  reference  species  is  obvious  onJ examination  of  the 
matrix  elements  in  the  text  in  the  discussion  of  the  matrix  manipulations 
required  for  transformation  of  reference  species.  As  stated  above,  such 
manipulations  do  not  affect  the 
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APPENDIX  II 

THE  INTRODUCTION  OP  REAL  GAS  EFFECTS 

The  expressions  derived  for  the  calculation  of  the  species  concentrations 
and  those  developed  for  the  calculation  of  the  properties  of  the  mixture  need 
to  he  modified  to  include  real  gas  effects.  In  principle,  these  modifications 
should  include  both  "internal"  effectB,  i.e.  modifications  of  the  internal 
energy  levels  of  individual  species  and  "external"  effects,  i.e.  those  due  to 
interactions  between  particles.  For  particle  densities  of  interest  here, 
however,  changes  in  the  internal  energy  levels  of  individual  species  will  be 
negligible.  As  a  result,  the  thermodynamic  properties  of  individual  species 
can  be  written  as  the  sum  of  an  ideal  gas  part  (calculated  as  indicated  in 
the  text)  and  a  part  due  to  particle  interactions  involving  the  species.  In 
this  appendix  we  shall  be  concerned  with  the  development  of  certain  correction 
terms  to  the  ideal  gas  expressions  developed  in  the  text. 

According  to  the  formalism  developed  for  the  ideal  gas,  the  inclusion  of 
expressions  for  the  species  activity  coefficients  will  automatically  extend 
the  calculations  of  the  species  concentrations  to  include  real  gas  effects. 

Such  expressions  need  to  be  included  in  equation  (10)  for  the  free  energy 
minimization  method  and  in  (17a)  for  the  equilibrium  constant  method.  The 
extension  of  the  calculation  of  the  thermodynamic  properties  of  the  mixture 
to  include  real  gases  is  then  completed  on  adding  to  the  relations  (23)  et  seq. 
developed  for  the  mixture  properties  of  an  ideal  gas  appropriate  terms  for 
real  gas  effects  and  using,  in  the  resulting  expressions,  concentrations  cal¬ 
culated  using  real  gas  activity  coefficients. 

I.  Real  Gas  Effects  on  the  Calculation  of  the  Species  Concentrations 

The  chemical  potential  of  a  species  in  a  gas  mixture  can  be  written 

U^T.P.x)  =  u°  (T)  +  InF  +  lnx1  +  E(AkV±)/RT 
RT  RT 


where  x  indicates  a  dependence  on  the  concentrations  of  all  species  and  where 
each  Ay  is  a  different  additive  effect.  An  activity  coefficient  can  be 
defined  for  each  of  these  effects  through  the  relation 


lny^  =  Aky/RT 


leading  to  an  expression  for  the  chemical  potential 
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y  , 
RT 


ui  W  +  InP  +  lnx  +  In  [Ify  ] 
RT  1  k  1 


The  expression  {10)  on  which  the  calculation  of  the  concentrations  in  the  free 
energy  minimization  method  is  based  then  need  only  be  modified  by  taking 
Ck) 

y  “  Ry,  including  as  many  of  the  k  effects  as  desired.  The  modification 
k 

of  (17a)  is  equally  straight  forward  merely  requiring  the  definition  of  a  y' 
for  each  of  the  k  effects.  It  should  be  noted  that  these  effects  can  include 
such  diverse  items  as  quantum  effects,  higher  virial  coefficients,  etc.  In 
short,  the  formalism  includes  all  effects  which  produce  additive  terms  In  the 
chemical  potential. 

The  effect  of  second  virial  coefficients  can  be  included  immediately  on 
writing  down  the  appropriate  correction  to  the  chemical  potential.  Thus 


so  that  InY 

i 


2  1  r 

RT'“  V  l  ^  kBik 
k=l 


GU> 


This  can  be  substituted  directly  into  (10)  for  the  free  energy  minimization 
method.  For  the  equilibrium  constant  method,  it  is  necessary  to  calculate 


[h  Yjij]/Y^ 


for  the  1th 


chemical  reaction. 


Thus 


y± 


exp 


1  v 


i 

1 

k- 


CkBik  J 


ao  that 


V  ^  2  1  2  * 

“  exP  ‘v  l  Ifk  Bik  "  V  B=1ck  V 
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and  Y 


exp[ 


p/c0  I 


V 


ij 


Z 

z  c 

k=l 


k 


l 

z  c 

k=l 


k 


Bik' 


CB2) 


where  V  .is  the  volume  of  one  mole  at  standard  conditions,  It  should  be 
noted  t$at  this  y'  should  be  substituted  into  (19)  as  a  factor  multiplying 

V 

/ 

The  limiting  law  Debye-Huckel’ correction  to  the  calculation  of  the 
species  concentrations  can  also  be  included  in  a  straightforward  manner.  In 
this  case 


where  N  is  the  number  of  particles  per  mole  at  standard  conditions,  V  the 
corresponding  volume,  D  the  dielectric  constant  of  the  mixture  and  Zg  ?he  ionic 
charge  of  each  ion  of  type  s.  It  follows  then  that 


lnYi 


Cb3) 


Substitution  of  this  expression  into  (10)  guarantees  the  inclusion  of  the 
limiting  law  Debye-Huckel  correction.  In  the  equilibrium  constant  method, 
this  correction  is  inserted  by  means  of  the  additional  factor 


r 

1 


IT  Y^ij 


exp 


(B4> 


Where  the  reactions  are  vjritten  in  terms  of  complete  ionization  to  ions  and 
electrons,  this  can  be  simplified  slightly  by  making  use  of  the  identity 
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l 

i 


v  Z2 

U  i 


Z. 

J 


for  that  case.  This  follows  from  the  fact  that  for  such  reactions  v 
except  for  the  electrons  while  for  the  electrons  v .  = 
and  Z,=l  the  electronic  valence. 


-Zj,  the  ionic 


ij 


=  0 

valence 


Woolley  has  developed  additional  corrections  for  various  other  effects 
including  third  virial  coefficients,  the  detailed  calculation  of  the  dielectric 
constant  and  the  finiteness  of  the  ionic  core. 


At  this  point,  the  advantages  of  the  direct  search  method  should  become 
apparent.  We  have  defined  direct  search  methods  to  include  all  those  for 
which  only  the  objective  function  needs  to  be  evaluated.  The  addition  of  y' 
factors  in  (17a)  as  complex  as  (B2)  and  (B4)  or  of  Inj  terms  in  (10)  of  the^ 
complexity  of  (Bl)  and  (B3)  would  require  much  more  complex  reprogramming  in 
methods  involving  the  use  of  derivatives  than  is  required  in  methods  for  which 
only  the  objective  function  is  evaluated.  Furthermore,  the  former  requires 
some  additional  mathematical  analysis  in  these  cases  whereas  the  latter 
approach  does  not. 

II.  Real  Gas  Effects  on  the  Calculation  of  the  Properties  of  the  Mixture 

In  part  I  of  this  appendix,  we  showed  how  the  species  concentrations  can 
be  calculated  including  various  real  gas  effects.  These  are  the  concentrations 
which  must  now  be  used  in  the  expressions  for  the  thermodynamic  properties  of 
the  mixture  in  terms  of  the  appropriate  properties  of  the  individual  species . 
These  expressions,  however,  are  not  simply  those  derived  in  the  text  for  the 
ideal  gas  mixture.  In  addition  to  these  ideal  gas  terms  there  must  be  included 
terms  which  specifically  refer  to  the  real  gas  effects  on  the  properties 
themselves.  Thus,  the  compressibility  factor  must  now  be  written 


Z  = 


PV 

RT 


Z*  +  Z  iZ 


where  Z*  is  the  ideal  gas  value,  i.e.  Z*  =  E  C.  and  the  &,  Z  are  expressions 

■  T  '  1  K 

1=1 

for  the  various  real  gas  effects.  Thus,  for  the  effect  of  the  second  virial 
coefficient 


hZ 


,  1  1 

„  -1  p/p  £  £  C  C  R 

V  °a=l  g-1  a  6 
o 


(B5) 
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where  B  „  is  the  second  virial  coefficient  which  describes  the  interaction 

ft 

between  species  oi  and  species  6.  For  the'  limiting  law  Debye-Huckel  effect 


AZ 


3 

£ 


(DkT)3/2 


(B6) 


where  Z 

a 


is  the  ionic  charge  of  species  a  . 


The  internal  energy  needs  now  to  be  computed  using 


E  e!  +  *  Ak  (  M 
RT  **  RT  k  \RT/ 


where  E*  is  the  ideal  gas  value  and,  for  the  effect  of  the  second  virial 
coefficient, 


A 


E_ 

RT 


A 

Z 

a=l 


C0T 


(B7) 


While  for  the  Debye-Huckel  effect,  the  correction  for  the  internal  energy 
is  three  times  (B6) ,  i.e.  three  times  that  for  the  compressibility  factor. 
The  entropy  is  to  be  calculated  from 

S  S*  +  Z  A./s  \ 

R-R-  k  k\&) 


where  S*  is  the  ideal  gas  value  where  the  second  virial  effect  ia  given  by 

r 


A 


J3 

R 


Z 

Z 

6=1 


C 

a 


Bag  +  T 

dT 


(B8) 


and  the  Debye-Huckel  limiting  law  effect  ^  given  by  (B6) . 

Corrections  are  easily  derived  for  the  various  properties  which  depend 
on  derivatives  of  these  properties  by  differentiation  of  (B5)  -  (B8)  as 
required. 
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